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SECTION  1 
INTRODUCTION 


This  report  describes  the  results  of  a  study  into  nonlinear 
system  identification  techniques.  The  primary  objective  was  to 
refine  previously  proposed  techniques  and  make  actual  circuit 
measurement  easier. 

The  need  for  nonlinear  system  identification  is  of  impor¬ 
tance  because  of  the  increasingly  dense  environment  of  present 
and  future  Air  Force  communication  systems.  Accurate  nonlinear 
system  identification  would  help  in  the  accurate  prediction  of 
such  nonlinear  effects  as,  for  example,  intermodulation  and  har¬ 
monic  distortion,  thus  leading  to  the  determination  of  RF  sus¬ 
ceptibility  of  various  electronic  equipment. 

A  promising  system  identification  technique  known  as  the 
pencil  of  functions  method  was  originally  proposed  by  Jain  [1974] 
and  later  applied  by  Ewen  [1975,1979]  to  the  identification  of 
the  nonlinear  transfer  functions  of  a  class  of  nonlinear  sys¬ 
tems.  Ewen's  identification  was  based  on  the  fact  that,  for 
lumped  parameter  circuits  containing  zero-memory  non  linear  i  ties 
between  circuit  nodes,  the  nonlinear  transfer  functions  have 
poles  which  are  uniquely  determined  by  the  poles  of  the  linear 
part  of  the  circuit.  Such  nonlinear  transfer  functions  can 
therefore  be  identified  by  finding  the  linear  transfer  function 
and  thus  its  poles  and  then  identifying  the  residues  or  zeros  of 
the  nonlinear  transfer  functions. 

Unfortunately,  the  identification  technique,  as  proposed  by 
Ewen  was  not  readily  implementable  in  practice  and  numerous  ques¬ 
tions  remained  as  to  the  appropriate  parameters  to  be  used  and 
the  performance  to  be  expected  in  a  realistic,  noisy  environment. 
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■  In  fact,  Ewen  found  that  for  a  two  pole  system,  his  re- 

H  quired  sampling  rate  was  4  to  10  times  the  highest  frequency  of 

■  the  passband  of  the  system.  The  sampling  rate  was  found  to  rise 
unreasonably  high  as  the  number  of  poles  was  increased.  Ewen 
concluded  that  current  A/D  device  availability  limited  the  band¬ 
width  of  two  pole  systems  that  could  be  analyzed  to  10  or,  at 

]  most,  30  kHz.  As  the  number  of  poles  increased,  the  sampling 

rate  requirements  increased.  Systems  with  more  than  two  poles 
required  a  sampling  rate  which  exceeded  the  Nyquist  sampling 
rate.  Thus,  additional  poles  limited  further  the  bandwidth  of 
;  the  systems  which  could  be  identified. 

Ewen  noted  that  the  fastest  commercially  available  A/D  con¬ 
verter  with  16  bit  resolution  had  the  capability  of  sampling  a 
signal  at  125  ksample/sec.  He  concluded  that,  at  least  at 

present,  only  two-pole  systems  with  a  maximum  bandwidth  of  up  to 
10K  or  30  kHz  could  be  analyzed.  For  a  4-pole  system,  the  reso¬ 
lution  required  in  an  A/D  converter  is  at  least  24  bit  per  sample 
at  125  ksample/sec  in  order  to  achieve  a  minimum  performance. 
This  again  indicates  that  at  the  present  time,  this  type  of  non- 
i  linear  system  analysis  would  be  limited  to  two-pole  systems. 

The  objective  of  our  study  was  to  solve  the  implementation 
problems  found  by  Ewen.  The  organization  of  this  report  is  as 
follows:  the  problem  is  first  defined  in  Chapter  2.  The  pencil- 

of-functions  method  and  Ewen's  technique  are  then  summarized  in 
Section  3.1.  The  identification  parameters  and  practical  diffi¬ 
culties  associated  with  Ewen's  method  are  analyzed  in  Section 
3.2.  Based  on  these  factors.  Section  3.2  concludes  with  an  out¬ 
line  of  the  detailed  objectives  of  the  present  effort. 

In  order  to  meet  these  objectives,  a  new  discrete  iterative 
identification  technique  is  defined  in  Chapter  4.  It  has  the  ad¬ 
vantage  of  being  defined  totally  on  discrete  samples  of  the  input 
and  output  processes,  thus  eliminating  a  major  difficulty  asso- 
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ciated  with  Ewen ' s  approach  which  stems  from  the  numerical  evalu¬ 
ation  of  repeated  continuous  integrals.  Moreover,  the  new  method 
permits  evaluation  of  the  reliability  of  the  identified  unknown 
poles  by  also  identifying  at  each  step  the  known  input  poles. 
The  distance  between  the  identified  input  poles  and  the  true  in¬ 
put  poles  is  indicative  of  the  accuracy  of  the  identified  unknown 
system  poles.  In  the  identification  of  wideband  systems,  itera¬ 
tion  proceeds  by  starting  with  the  identification  of  the  lowest 
frequency  poles,  then  using  previously  identified  lower  frequency 
poles  as  knowns  to  identify  the  next  higher  frequency  poles. 

The  problem  of  identifying  a  system  in  the  presence  of 
noise  is  examined  in  Chapter  5.  The  effects  of  noise  are  inves¬ 
tigated  and  desirable  identification  parameters  are  specified  us¬ 
ing  simulation  results. 

Results  of  simulations  of  the  discrete  iterative  identifi¬ 
cation  procedure  are  presented  in  Chapter  6  for  wideband  and  nar¬ 
rowband  systems  containing  up  to  four  closely  grouped  moles.  The 
simulations  were  performed  with  varying  levels  of  noise  added  to 
the  system  output.  The  standard  deviation  of  the  noise  ranged 
from  0  to  5%  of  the  maximum  system  output.  The  highest  sampling 
frequency  used  is  not  greater  than  the  Nyquist  rate.  In  all 
cases  the  correct  number  of  poles  was  identified.  The  average 
error  in  pole  location  ranged  from  .9%  for  the  noiseless  case  to 
15%  for  additive  noise  with  a  standard  deviation  equal  to  5%  of 
the  maximum  value  of  the  system  output.  The  feasibility  of 
achieving  identification  of  this  accuracy  at  sampling  rates  no 
higher  than  the  Nyquist  rate,  greatly  simplifies  future  implemen¬ 
tation.  The  new  method  substantially  decreases  the  required 
sampling  rate  and,  because  of  the  greater  tolerance  to  noise,  the 
necessary  A/D  converter  resolution. 

Possible  performance  assessement  criteria  are  outlined  in 
Chapter  7.  Two  types  of  criteria  are  possible.  In  Section  7.2, 
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the  accuracy  of  the  identified  poles  and  zeros/residues  are  as¬ 
sessed  directly.  In  Section  7.3,  global  measures  of  the  accuracy 
of  the  identified  impulse  response  or  transfer  function  are  dis¬ 
cussed.  It  should  be  noted  that  nonlinear  transfer  function 
identification  makes  use  of  the  results  of  the  linear  system 
identification,  but  requires  the  additional  identification  of  the 
zeros/residues  of  the  nonlinear  transfer  functions  which  was  not 
undertaken  during  the  current  effort.  Thus,  the  impact  of  the 
linear  system  parameter  and  global  error  measures  on  the  quality 
of  nonlinear  transfer  function  identification  was  not  evaluated 
in  this  study. 

Finally,  Chapter  8  summarizes  the  characteristics  of  the 
discrete  iterative  identification  technique  and  presents  conclu¬ 
sions  as  to  its  expected  performance.  The  main  characteristics 
of  the  postulated  technique  are: 

•  Completely  discrete  formulation. 

•  Selective  use  of  the  knowledge  of  the  input  to 
test  for  regions  where  poles  are  present  and  to 
determine  the  reliability  of  the  identified 
poles. 

•  Iterative  decrease  of  the  sampling  interval  T. 

•  Use  of  previously  identified  poles  as  knowns  in 
the  set  defining  the  Gram  determinant  during  the 
identification  of  additional  poles. 

The  advantages  of  the  discrete  iterative  approach  can  be 
summarized  as  follows: 

•  It  does  not  require  repeated  integrations;  in¬ 
stead,  time  shifts  are  used. 


It  permits  the  determination  of  regions  where 
system  poles  are  present. 


•  It  aids  in  determining  the  number  of  poles  being 
iden  tif ied. 


•  The  accuracy  of  an  identified  input  pole  when  it 
is  far  from  system  poles  gives  a  measure  of  the 
best  accuracy  that  can  be  expected. 

•  Identifying  the  input  poles  while  using  the  iden¬ 
tified  poles  as  knowns  gives  a  measure  of  the  ac¬ 
curacy  of  the  identified  poles. 

•  The  use  of  identified  poles  as  knowns  is  espe¬ 
cially  useful  for  wideband  systems. 

•  Linear  identification  requires  the  measurement  of 
only  approximately  twenty  output  samples. 

•  The  required  sampling  frequency  is  not  higher 
than  the  Nyquist  rate. 


Based  on  simulation  results,  it  is  conservatively  concluded 
that  the  method  should  be  capable  in  practice  of  identifying  sys¬ 
tems  containing  separate  groups  of  four  or  fewer  closely  packed 
poles  in  the  presence  of  additive  output  noise  whose  standard  de¬ 
viation  is  of  the  order  of  1%  of  the  maximum  system  signal  out¬ 
put.  If  the  output  noise  is  mainly  due  to  quantization,  this 
should  permit  the  use  of  9  or  10  bit  quantizers.  The  required 
sampling  frequency  is  approximately  equal  to  the  frequency  of  the 
highest  system  pole.  These  values  should  permit  a  relatively 
straightforward  implementation. 


SECTION  2 

PROBLEM  DEFINITION 


In  the  increasingly  dense  communications  environment  found 
in  many  defense  communication  systems,  it  is  of  great  interest  to 
be  able  to  predict  the  amount  of  intermodulation  and  harmonic 
distortion  generated  in  these  various  equipments.  This  entails 
accurate  descriptions  of  the  nonlinear  characteristics  of  such 
equipment. 

In  the  case  where  the  nonlinear  behavior  exhibits  no  jumps 
or  hysterisis,  the  nonlinear  system  can  be  represented  by  a  Vol- 
terra  series.  Its  output  y(t)  due  to  excitation  x(t)  is  then 
given  by 


CD  CO  CO  pj 

y ( t )  =  l  /  ...  /  h  (u  ,...,u  )  n  x(t-u.)du. 

n  =  l  -“-<■»  i  =1 


«  l 


n=l 


yn(t) 


(2.1) 


where,  yn(t)  is  nth-order  system  output.  The  system  is 

then  characterized,  that  is  its  input-output  relation  is  com¬ 
pletely  specified  by  the  nth-order  impulse  responses 
hn(u1,...  ,un),  n=l , 2 , . . .  . 

If  the  system's  nonlinearity  is  mild,  the  output  is  given 
by  the  first  few,  predominantly  three  first  terms  of  the  series: 

3 

y(t)  =  l  y  ( t ) .  (2.2) 

n=l  n 
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In  order  to  characterize  such  a  system,  it  is  necessary  to  deter¬ 
mine  hj^u^),  h2<u^,U2),  h^  (u  ,  U2 ,  Uj )  or-  equivalently,  the 
higher-order  transfer  functions 


H1(s)  =  £1[h1(u)] 

h2(si's2)  =  £2fh2^Ul'U2^ 
s 1 ' s 2 ' S3^  =  £ 3^ ^3 ^ u 1 ' u 2 ’ U 3 ^ 


(2.3) 


where  £n  denotes  n-dimensional  Laplace  transformation. 

This  problem  is,  in  general,  extremely  complex  since  it  in¬ 
volves  the  determination  of  multidimensional  functions.  It  has, 
however,  been  shown  that,  in  the  case  where  the  nonlinear  system 
is  a  lumped  parameter  circuit  with  zero-memory  nonlinearities  be¬ 
tween  circuit  nodes,  the  equivalent  2nd  and  3rd-order  transfer 
function  poles  can  be  obtained  from  the  poles  of  the  linear 
transfer  function  [Graham  and  Ehrman  (1973);  Ewen  (1975)]. 


More  precisely,  let  the  transfer  function  of  the  linear 
part  of  the  system  be  given  by 


H^s) 


M 


n  (s+q.) 
i=l  1 


N 


M<N 


n  (s+p.) 
i=l  1 


N  R. 

=  l  1 


(2.4) 


i=l  S+P] 


where ,  — q  ,  ■*  q  “Pi  •  ~  P2/«»./— ?  ^1 ' ^2  f  •••*  ,  are  , 
respectively,  the  zeros,  poles  and  residues  of  the  transfer  func¬ 
tion  and  it  has  been  assumed  for  notational  simplicity  that  the 
poles  are  distinct. 


Then,  the  2nd  and  3rd-order  transfer  functions  are  given, 
respectively,  by  [Ewen  (1975)] 


H2^  sr  s2^ 


L  N 

l  I 

k1=l  k  =1  k1k2 


Sl+S2+2ek2 

[  s1+s2+aki+ak2  j  (  s2+ak2)  (  s^l 


(2.5) 


H3^s1's2'S3^  ~ 


,  J  L  N 

4  I  l  l 


klk2k3 


k1=l  k2=l  k3=l  ls1+s2+s3  +  (ak1+ak2+ak3)^ 


, _  1 _ ^  + _ 1  _ r 

^sl+ak3  ^sl+s3+^ak2"ak3  ^  )  ^s2+ak3^sl'*'s2+^ak2+ak3  ^ 


+  1  ,  _  1  _ 

(  Sl  +ak^)  (  s1+s2  +(  aR ^+ak  )  (  s3+ak^)  (  s2+s3+(  ak  +aj^TT 


1  _  +  _  1  _ 

(  s2+ak3^  (  s2+S3+^ak2+ak3  J  J  ^S3+ak3^Sl+S3+^ak2+ak3^  J 


where,  the  quantities  J,  L,  a^ ,  are  uniquely  determined 
by  the  poles  of  the  linear  transfer  function  H^fs).  The  com¬ 
plete  identification  of  the  nonlinear  transfer  functions  of  the 
system  requires,  therefore,  the  identification  of  the  linear 
transfer  function  and  the  determination  of  the  constants 

Ak  k  <  Ck  k  k  t*le  permissible  values  of  k^,  k£»  and  kj. 

12  12  3 

The  poles  of  the  linear  transfer  function  thus  play  a  cru¬ 
cial  role  in  the  identification  of  the  linear  and  nonlinear  parts 
of  the  system.  They  not  only  specify  the  denominator  of  H^s) 
but  linear  combinations  of  these  poles  determine  the  poles  of 
f^fs^  ,s2  )  a  nd  H^(S2/S2,S2)» 
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In  the  next  section,  an  identification  technique,  known  as 
the  pencil  of  functions  method,  developed  by  Jain  [Jain  (1974)] 
and  applied  by  Ewen  [Ewen  (1975),  (1979)]  to  the  above  problem  is 
presented.  The  critical  parameters  of  the  technique  and  the 
practical  difficulties  involved  in  applying  it  are  isolated  and 
analyzed.  An  efficient  and  practical  method  for  the  identifica¬ 
tion  of  the  poles  of  a  linear  system  is  then  developed  and  illus¬ 
trated  by  simulation  examples. 


SECTION  3 

PENCIL  OF  FUNCTIONS  METHOD 


3 . 0  Summary  of  the  Method 

A  promising  technique  for  determining,  from  input/output 
measurements,  the  poles  of  the  transfer  function  of  a  linear  sys¬ 
tem  has  been  developed  by  Jain  {1974],  Known  as  the  pencil  of 
functions  method,  this  technique  was  then  applied,  by  Ewen  [1975, 
1979]  to  the  identification  of  the  2nd  and  3rd-order  transfer 
functions  of  a  Volterra  system.  We  next  outline  the  pencil  of 
functions  method  as  originally  presented  by  Jain  and  used  by 
Ewen.  This  technique  was  the  starting  point  of  the  present  con¬ 
tract.  Subsequent  developments  by  Jain  and  Osman  [1979]  under 
separate  contract  No.  F30602-75-C-011 8,  conducted  in  parallel 
with  the  present  contract,  are  presented  where  appropriate. 

Suppose  that  the  transfer  function  of  the  linear  part  of 
the  system  is  given  by 


H  (  s  )  = 


n  (s+q.) 
i=l  1 
N 

n  (s+p.) 
i=l  1 


i=l  s+pi 


_  Q(s) 

P(s) 


(3.1) 


Excite  the  system  by  x^(t)  .  With  x^(t)  known  or  measured, 
measure  the  system  output  y^t).  By  successive  integrations, 
compute 


Ewen  selected  x^(t)  exponential  for  ease  of  analysis 
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(3.2) 


t 

x .  (t )  =  /  x .  . (u  )du 

1  0 

t 

y.  (t)  *  /  y,  ,  (u)du,  i  =  l,2  . 

1  0 


Suppose  that  N,  the  number  of  poles,  is  known.  Form  the  set 

(3.3) 

i(y1(t)-^y2<t))  »(y2(fc>-*y3(t))  »•  •  •  ^ ' 


x2(t ) , .  .  .  *xN+1  (t  )} 

where  X  is  constant.  Test  the  set  for  linear  dependence  by 
forming  the  linear  combination 


c1(y1(t)-Xy2(t))  +  ...  +  cn(yN<t)-XyN+1(t)) 
+  d1x2(t)  +  ...  +  dNxN+1<t)  =  0 


(3.4) 


The  set  is  linearly  independent  if  the  only  solution  to  (3.4)  is 

ci  =  c2  =  ...  =  cN  =  d^  =  ...  =  dN  =  0. 

It  is  linearly  dependent  otherwise. 

Taking  the  Laplace  transform  of  (3.3)  an  equivalent  test 
for  linear  independence  is 


(3.5) 

C1(Y1(S)  -  X Y 2 ( s ) )  +  ...  +  Cn(YN(s)  -  XYN  +  1(S)) 

+  d |X  2 ( s )  +  ...  +  dNXN+1(s)  =  0 
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But,  using  (3.1)  and  (3.2), 


Ylls)  -  Hffxl(s) 

X.(s)  =  X1(s)  /  s1”1  (3.6) 

Y^s)  =  Y 1  (  s  )  /  s1"1,  i=l,  2,  ...  . 

The  test  for  linear  independence  in  (3.5)  therefore  becomes 

Q(sl  (S-UfCjS*1  1  +  CjS  2*  •  ••  *  cN)  7) 

+  P(s)(dlsN'1+  ...  +  d„)  =0 

The  left  hand  side  of  (3.7)  contains  2N  coefficients.  If 
(s-X)  is  not  a  factor  of  P(s)  that  is,  if  X  is  not  a  pole  of 
H  ( s )  the  polynomial  in  (3.7)  is  of  degree  2N-1,  containing 
terms  in  1,  s,  ...,  s2N-^,  and  it  has  2^  coefficients.  The 
set  is  therefore  linearly  independent. 

If  x  is  a  pole  of  H(s)  or,  equivalently,  (s-X)  is  a 
factor  of  P(s), 

P(s)  =  p' (s) (s-x)  (3.8 ) 

where,  p'(s)  is  of  degree  N-l.  Dividing  (3.7)  by  (s-X)  '  re¬ 
sults  in  a  polynomial  of  degree  2N-2.  Since  there  are  still 
2N  coefficients,  the  set  is  linearly  dependent. 

It  follows  that  the  poles  of  H(s)  are  the  values  of  X 
for  which  the  set  in  (3.3)  is  linearly  dependent.  The  linear  de¬ 
pendence  of  a  set  of  time  functions  can  be  readily  checked  using 
the  Gram  determinant  which,  for  the  set  of  functions 
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{to  ^  ( t )  ,  to  2 (  t )  ,  •  •  •  i  <o jc<  t )  } ,  0  <  t  <T 


(3.9  ) 


is  defined  as 


Gk  =  det 


W  j  •  •  • 


“lwk 


to  ^  •  •  •  (0  2^ 


(Okco1  .  .  .  (0kt0k 


(3.10; 


where  the  inner  product  (o^uk  over  the  interval  [0,T]  is  equal  to 


(0.(0-  =  /  u).(x)to.  (t)dT,  i  ,  j=l,2  ,  .  .  .  ,k 
A  J  n  A  J 


and  where  *  denotes  complex  conjugation. 

Hence,  the  poles  of  H(s)  are  the  values  of  X  for  which 
the  Gram  determinant  of  the  set  of  functions  in  (3.3)  is  equal  to 
zero. 

Equivalently,  it  can  be  shown  (Jain  (1974)]  that  the  poles 
are  the  roots  of  the  polynomial  in  X 


N 


l 

i=0 


x 


N-i 


^G2N+1  \  i+1,  i+1)^ 


1/2 


=  0 


(3.11) 


where,  [G2n+i1  ( i.+ 1  i+i)  *s  the  value  of  the  (i+l,i+l)th  co¬ 
factor  of  the  Gram  determinant  of  the  set 


{y]L(t),  ...  ,  yN+1(t),  x2(t). 


xN+1(t)}  .  (3.12) 
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In  the  above,  it  was  assumed  that  N,  the  number  of  poles, 
was  known.  In  practice,  this  number  must,  of  course,  be  found 
before  determining  the  poles  using  (3.11).  The  number  N  can  be 
determined  by  considering  the  set 

l y 1 ( t ) ,  ...,  y  ( t ) ,  x  2  ( t ) ,  • • • ,  x  ^  ( t )  ]  .  (3.13) 

Using  the  same  approach  as  above,  this  set  is  linearly  independ¬ 
ent  if  K<N;  it  is  linearly  dependent  if  K>N+1.  The  number  of 
poles  can,  therefore,  be  found  by  calculating  the  Gram  determin¬ 
ant  for  the  set  in  (3.13)  for  increasing  values  of  K, 

I 

K=2 , 3 , . . .  .  If  K  is  the  smallest  K  for  which  the  Gram  de¬ 
terminant  is  zero,  that  is,  the  first  value  for  which  the  set  is 
linearly  dependent,  the  number  of  poles  is  equal  to 

N  =  k'  -  1  .  (3.14 ) 

3. 2  Identification  Parameters  and  Practical  Difficulties  Asso¬ 
ciated  with  Ewen's  Methods 

The  identification  of  an  unknown  circuit  using  the  pencil 
of  functions  approach  as  proposed  by  Ewen  [Ewen  (1975),  (1979)] 

and  as  described  in  Section  3.1  therefore  involves  the  following 
steps: 

(1)  Excite  the  system 

(2)  Measure  the  system  outputs  and  inputs 

(3)  Calculate  X2(t),...,  xk(t),  y2<t),  ...,  y^ft)  as  de¬ 

fined  in  (3.2)  by  successive  integrations. 


(4)  Calculate  the  inner  products  over  an  interval  T. 

(5)  Calculate  the  Gram  determinant  of  the  set  of  functions 
in  step  (3)  for  increasing  values  of  K. 

(6)  Determine  the  number  of  poles  N  as  K'-l,  where  K'  is 
the  smallest  value  of  K  for  which  the  Gram  determinant 
is  zero. 

(7)  Find  the  poles  using  (3.11). 

Because  of  the  amount  of  computation  which  is  necessary,  it 
is  more  efficient  to  perform  the  calculations  using  a  computer. 
This  suggests  the  iden tif ication  set-up  shown  in  Figure  3.1  where 
the  double  arrow  at  the  input  indicates  that  the  input  can  either 
be  generated  in  analog  form  and  then  converted  to  digital  form  or 
generated  by  the  computer  and  converted  to  analog  form. 

Note  that,  in  addition  to  the  blocks  shown  in  Figure  3.1,  a 
practical  identification  implementation  would  necessitate  use  of 
a  linear  output  and  input  amplifier /attenuator  with  large  dynamic 
range  to  insure  that  the  system  input  is  of  the  proper  level  and 
that  y(t)  uses  the  full  range  of  the  analog  to  digital  con¬ 
verter.  The  contract  under  which  the  present  effort  was  per¬ 
formed  called  originally  for  the  design  and  the  evaluation  of  a 
practical  implementation  of  the  above  technique.  This  was  how¬ 
ever  abandoned  because  of  the  following  practical  difficulties 
and  unresolved  questions  associated  with  the  method: 

(1)  What  is  the  impact  of  using  different  excitations?  In 
the  original  technique  only  a  decaying  exponential  in¬ 
put  was  considered. 
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T  ■  * 


(2)  Because  of  the  digital  implementation  the  output  sig¬ 
nal  must  be  sampled  and  quantized.  A  sampling  rate 
must  therefore  be  specified  which  will  assure  optimum 
performance  and  which  will  be  practical  with  existing 
analog  to  digital  converters. 

(3)  Simpson's  rule  was  used  in  evaluating  the  integration 
of  a  function  as  given  by 

i_  _ 

[y(0)  +  4y  ( A  t )  +  2y(2At)  +  4y(3At)  + 

...  +  2y((2n-2)At)  4y((2n-l)At)  +  y  ( 2nAt )  ]  . 


b 

/  y(t )dt  = 
a 


Each  integration  reduces  the  number  of  samples  by  one 
half  and  therefore  increases  the  error.  As  the  number 
of  poles  of  the  system  increases,  the  number  of  re¬ 
quired  integrations  increases  accordingly.  Integra¬ 
tion  by  Simpson's  rule  introduces  error  at  each 
step.  After  a  number  of  integrations,  therefore,  the 
error  becomes  excessively  large.  The  accuracy  of  suc¬ 
cessive  integrals  decreases  rapidly.  This  loss  of  ac¬ 
curacy  due  to  numerical  integration  can  be  overcome  by 
increasing  bits  in  the  sample  and  the  sampling  rate. 
However,  the  specifications  of  hardware  devices  which 
implement  analog  to  digital  (A/D)  and  digital  to  ana¬ 
log  (D/a)  conversions  and  the  high  rate  of  sampling 
rapidly  exceed  the  range  of  current  technology  as  the 
number  of  poles  of  the  system  increases.  For  example, 
Ewen  states  that  if  we  were  to  analyze  a  system  of 
four  poles,  the  A/D  converter  would  be  required  to 
have  20-24  bits  of  resolution.  The  present  commer¬ 
cially  available  A/D  converters  are  capable  of  at  most 
a  16  bit  resolution. 
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The  performance  of  successive  numerical  integra¬ 
tions  has  a  large  impact  on  the  necessary  sampling 
rate.  Ewen  concluded  that  the  sampling  rate  required 
is  determined  by  the  need  for  accuracy  of  numerical 
integration  and  not  by  Nyquist  sampling  constraints. 
For  a  two  pole  system,  the  required  sampling  rate  is  4 
to  10  times  the  highest  frequency  of  the  passband  of 
the  system  if  we  desire  a  mean  square  error  between 
the  actual  and  predicted  system  output  to  be  smaller 
than  the  mean  square  error  induced  by  a  10  percent 
error  in  each  pole  and  residue.  The  sampling  rate  may 
rise  unreasonably  high  as  the  number  of  poles  in¬ 
creases.  Considering  errors  in  identified  poles,  Ewen 
concluded  that  the  bandwidth  of  a  two  pole  system  that 
can  be  analyzed  is  limited  to  10  to  30  kHz  by  current 
device  availability. 

Ewen  noted  that  the  fastest  commercially  avail¬ 
able  A/D  converter  with  16  bit  resolution  had  the  cap¬ 
ability  of  sampling  a  signal  at  124  kHz.  This  re¬ 
sulted  in  his  conclusion  that,  at  least  at  present, 
only  systems  with  a  maximum  bandwidth  of  up  to  10  kHz 
or  30  kHz  could  be  analyzed. 

(4)  .The  inner  product  interval  T  which,  in  noisy  practical 

situations,  will  result  in  the  specified  error  in  the 
location  of  the  identified  poles,  cannot  be  readily 
determined. 

(5)  The  measured  digital  output  which  is  used  will  contain 
noise  due  essentially  to  three  sources:  quantization 
noise  from  the  A/D  converter,  thermal  noise  originat¬ 
ing  in  the  unknown  circuit  and  measurement  errors. 
The  inner  products  which  are  the  elements  of  the  Gram 
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matrix  will  therefore  be  noisy.  As  a  consequence,  the 
Gram  determinant  will  not  be  zero  for  any  k.  This 
will,  of  course,  make  it  difficult  to  determine  the 
number  of  system  poles.  Even  if  the  correct  value 
of  N  is  determined  errors  in  the  inner  products  will 
result  in  errors  in  pole  locations.  These  errors  in 
pole  locations  can,  of  course,  be  expected  to  be  in¬ 
creased  if  the  wrong  value  of  N  is  chosen. 

(6)  Quality  criteria  must  be  studied  which  will  permit  to 
quantify  the  performance  of  the  identification  tech¬ 
nique  and,  after  the  global  identification  of  the 
system  (the  idenfci f ication  of  the  poles  arc 
zeros/residues,  traded  off  against  the  problems 
identified  in  (l)-(5). 

After  the  initial  study  resulted  in  the  above  conclusions, 
it  was  decided  with  the  concurrence  of  RADC  to  redirect  the  ef¬ 
fort  to  the  study  of  the  aforementioned  difficulties  and  unre¬ 
solved  questions.  In  Section  4,  we  describe  the  identification 
procedure  that  was  developed  to  combat  these  problems.  More  pre¬ 
cisely,  the  following  items  were  addressed: 


Selection  of  excitation  waveforms  and  techniques. 

Selection  of  input  and  output  functionals  to  be 
used  in  forming  the  Gram  determinant. 

Determination  of  correlation  intervals  to  be  used 
in  the  calculation  of  the  inner  products. 

Investigation  of  stopping  rules  for  the  determin¬ 
ation  of  the  number  of  linear  system  poles  in  the 
presence  of  system  noise,  quantization  noise  and 
as  a  function  of  the  relative  importance  of  the 
poles. 
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Evaluation  of  the  effectiveness  of  the  proposed 
techniques. 

Analysis  of  the  expected  errors  in  pole  locations 
due  to  quantization  noise,  underspecification  of 
the  number  of  poles  and  overspecification  of  the 
number  of  poles. 

Assessment  of  the  class  of  systems  that  can  be 
identified  given  practical  sampling  rates  and 
quantization  accuracies. 

Study  of  measurement  and  computational  techniques 
that  minimize  errors  in  pole  locations  by  averag¬ 
ing  out  results  of  several  measurements  or  of 
several  stages  of  the  identification  process. 
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SECTION  4 

DISCRETE  ITERATIVE  APPROACH 


4 . 1  Principle  of  the  Discrete  Method 

4.1.1  Basic  Method 

As  discussed  in  Section  3.2,  because  of  the  noise  intro¬ 
duced  by  the  numerical  calculation  of  repeated  integrals  and  the 
ensuing  necessity  of  using  very  high  sampling  rates,  a  primary 
objective  of  any  practical  implementation  is  to  modify  the  pro¬ 
cedure  so  as  to  eliminate  this  difficulty. 

This  can  be  accomplished  by  analyzing  the  problem  directly 
in  the  discrete  sampled  domain  instead  of  implementing  a  digital 
technique  which  emulates  analog  operations,  i.e.,  the  numerical 
computation  of  successive  integrals. 

As  in  (3.1)  suppose  that  the  linear  part  of  the  system  is 
characterized  by  input-output  transfer  function 

M 

n  (s+q.) 

H  ( s )  =  iip -  ,  M<N  (4.1) 

n  (s+p.) 
i=l  1 

N  R. 

r  l 


where,  for  notational  simplicity,  it  is  assumed  that  the  system 
poles  are  distinct.  If  the  unknown  circuit  is  excited  by  the  sum 
of  I  exponentials 
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X(t) 


i=l 


a  .  e 
1 


-bit 


(4.2) 


the  system  output  is  equal  to 


N  -p.t  I  -b.t 

y(t)  =  l  c.e  1  +  l  d.e  1 

i=l  i=  1 


(4.3) 


where,  the  constants  and  are,  respectively,  given  by 

I  R, 


Ci  A,  b.-p, 


j=l  “j 


N  R. 

di  =  A.  p^ET 


Note  that  the  assumption  that  the  input  is  the  sum  of  ex¬ 
ponentials  is  equivalent  to  the  assumption  that  the  Laplace 
transform  of  x(t)  is  the  ratio  of  two  polynomials,  in  s  with 
the  denominator  polynomial  of  higher  degree  than  the  numerator 
polynomial  and  the  denominator  polynomial  having  distinct 
roots.  Sampling  y(t)  with  sampling  frequency  fs=l/T,  sampl¬ 
ing  interval  T,  the  sampled  version  of  the  output,  over  an  in¬ 
terval  i 0, LT] ,  is  given  by 


yx(nT) 


N  -np.T  I  -nb-T 

£ c.e  +  £  d  .e  ,  0<n*L  . 

i=l  1  i=l  1 


(4.4  ) 


Note  that  the  choice  of  L,  the  number  of  samples  used,  and  of 
the  sampling  period  T  are  discussed  in  Section  5.3.2. 
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By  successive  shifts 


form 


y2(nT)  =  y^fn+DT) 


”  ~(n+l)p,T 

l  c.e  1 

i=l  1 


+ 


i  -(n+l)b.T 
l  d  e  i 

i*l  1 


N 


=  1 
i  =  1 


-n  p  .  T 
c,(2)e  1 


I 

I 

i  =  l 


-nb.  T 
di(2)e  1  , 


0<n<L 


yK(nT>  -  yK.1((n*1)T)  ,  yi((„+K-l,T) 


L  c.(k)e 
i=l  x 


-nPiT 


l  .  -np  T 
,1  d  (k)e  1  , 
i=l  1 


0<n<L 


where 


ci(k)  =  ci(k-l)e  PiT 

di<k>  -  di(k-l)e’biT 


c.e 


-(k-lJPjT 


~ (k-l)b.T 

d.e  1  V-o  o 

l  '  k-2,3, .  . .  ,K  . 
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Note  that  each  of  the  y^(nT),  i=l,2,...,K  is  a  linear  combina- 

-np.T  — nPj.T  -nb.T 

tion  of  the  functions  e  e  ,  e  1  , 

-nbjT 

e  .  Since  these  functions  are  linearly  independent  they 

can  be  considered  to  form  a  basis  for  y^(nT),  i=l,2,...,K. 

Consider  the  set  of  shifted  versions  of  the  sampled  output 


{y^(nT ) ,y2(nT) , . . . ,yK(nT)} 


(4.6) 


and  test  it  for  linear  independence  by  forming  the  linear  combi¬ 
nation 


alyl(nT)  +  a2y2(nT)  +  ***  +  aKyK(nT)  =  °'0<n<L  •  (4.7) 

-np  T  -npNT 

Since  the  linear  combination  in  (4.7)  has  e  ,  ...,  e  , 

-nb.T  nb^T 

e  1  ,  . . . f  e  as  its  basis,  it  is  a  (n+I )-dimensional 

function.  The  left-hand  side  of  (4.7)  contains  K  coeffi¬ 
cients.  It  follows  that  the  set  is  linearly  independent  if 

K<N+I  (4.8) 

and  that  it  is  linearly  dependent  if 

K>N+I  .  (4.9) 
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As  in  Section  3.1,  the  linear  dependence  can  be  checked  by 
using  the  Gram  determinant. 

The  set  of  functions  in  (4.6)  is  linearly  dependent  if  the 
Gram  determinant 


G^  =  det 


*2*1 


*1*2 


*2*2 


*'1*K 


*2*K 


*K* 1  -  *K*2 


*K*K 


(4.10) 


*  0 


where , 


_  ^ 

yHy^  =  l  y. (nT)y. (nT),  .i, j=l,2,. .. ,K  . 
13  n=0  1  3 


The  set  is  linearly  independent  if 
Gk  /  0. 

The  dimensionality  N+I  of  the  basis,  and  thus  the  number  N  of 
poles,  can  be  determined  by  calculating  the  Gram  determinant  for 
increasing  values  of  K.  If  K'  is  the  smallest  value  of  K 
such  that  Gk,=0,  then 
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N+I 


K  ’  - 1  . 


(4.11) 


Having  determined  N+I,  the  system  poles  -p^,  . ..,  -pN  and  the 
input  poles  -b^,  •••»  -bj  can  be  found  using  a  similar  ap¬ 
proach. 

Consider  the  set 


(4.12) 

{(Xy1(nT)-y2(nT)) ,(Xy2(nT)-y3(nT)) , . . . ,(XyN+J (nT)-yN+J+1 (nT) ) }  . 


Using  (4.5),  the  test  for  linear  independence  can  be  written  as 


N  -p.T  -np.T  I  -b.T  -nb.T 

l  c.(x-e  x)e  +  l  d.(\-e  )e 

i  i  A 


i=l 


i  =  l 


(4.13) 


+  a. 


N  ~P,T  -p.T  -np.T  I  -b.T  -b.T  -nb.T 
l  c.e  (x-e  ) e  +  I  d.e  (x-e  1  )e  1 

i=l  1  i=l  1 


+  •  •  •  +  a 


N+I 


?  cie-«M+I“1>PiT(X-e-PiT)e"nPiT 
i  =  l  1 


I  - ( N+I-l ) b . T  -b.T  -nb.T 

l  die  1  ( ^~e  1  ) e  1 


i  =  l 


=  0 


If  inX/T  is  not  one  of  the  system  or  input  poles,  that  is,  if 

-p.T 

X  *  e  1 

and 


i=l,2,. . . , N 


(4.14) 


-b.T 

X  *  e  ,  i=l , 2 , . . . , I  , 
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the  left-hand  side  of  (4.13)  is  of  dimension  N+I  and  contains 
N  +  I  coefficients.  It  follows  that  the  set  in  (4.12)  is  then 
linearly  independent. 

If  AnX/T  is  one  of  the  system  or  input  poles, 


-p.T 

X  =  e  , 

1<  i<N 

(4.15) 

or 

-b  T 

X  =  e  , 

1< i< I  , 

the  dimension  of  the  space  in  (4.13)  is  N+I-l  but  the  left-hand 
side  still  contains  N+I  coefficients.  It  follows  that  the  set 
is  linearly  dependent. 

Hence,  the  system  and  input  poles  can  be  determined  from 
the  values  of  X  for  which  the  set  in  (4.12)  is  linearly  depen¬ 
dent,  or  equivalently,  for  which  the  Gram  determinant  of  the  set 
if  zero.  Using  an  approach  due  to  Jain  11974]  as  in  (3.11)  it 
can  be  shown  that  these  values  of  x  are  the  roots  of  the 
polynomial 


N+I 


l 

i=0 


XN+!-! 


([GN+I  +  1-I  (  i  +  1, i+1 ) ) 


1/2 


0 


(4.16) 


where,  [G„. T . , 1 ,  ■  ,  ,  s is  the  (i+l,i+l)th  cofactor  of  the 
L  N+I+l J ( 1+ 1 r 1+ 1 ) 

Gram  determinant  of  the  set 


(yx(nT),  y2(nT),  •  ••,  yN+I+1CnTW  • 


(4.17) 
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If  ' ^2 ' *  *  * ' ^N+I  represent  the  N+I  roots  of  the  polynomial 
in  (4.16),  the  system  and  input  poles  can  be  gained  from 


AnXi/T,  i=l, 2, . . . ,N+I 


(4.18) 


Our  discrete  sampling  and  shift  identification  technique  is 
very  similar  to  the  pencil  of  functions  method  described  in  Chap¬ 
ter  3.  The  essential  differences  are  that 


•  The  proposed  technique  is  purely  discrete:  the 

successive  integrals  have  been  replaced  by  successive 
shifts. 

•  If  the  system  is  excited  by  the  sum  of  I  exponential 
excitations,  (N+I)  poles  are  identified  based  on  the 
determinant  of  a  (N+I +1 ) x (N+I+l )  matrix.  The  tech¬ 
nique  in  Chapter  3  used  a  ( 2N+1 ) x ( 2N+1 )  matrix  to 
identify  N  poles. 


Note  that,  in  the  above,  the  input  poles  are  identified  to¬ 
gether  with  the  system  poles.  In  the  following  section  the  iden¬ 
tification  method  derived  here  is  slightly  modified  to  make  use 
cf  the  knowledge  of  the  input  poles. 


4.1.2  Use  of  Input 

During  the  identification  of  an  unknown  circuit,  the  input 
to  the  circuit  is,  of  course,  known.  This  knowledge  can  there¬ 
fore  be  used  in  the  identification  technique. 

-nb.T  -nb.T  -nbjT 

Suppose  that  e  ,  e  ,  . .  .  ,  e  are  known 

and  define  the  set 
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-nb,T  -nbTT 

{ y ^ (nT ) ,y 2 (nT ) , • • • ,y^(nT ) , e  ,***,e  1  }  .  (4.19) 


Test  this  set  for  linear  independence  by  considering 

-nb,T 

alyl(nT)  +  a2^2(nT)  +  ***  +  aKyK^nT^  +  aK+le  (4.20) 

-nb  T 

+  •••  +  aK+^e  =  0  ,  0<n<L 


where  y^(nT),  i=l,2,...,K  is  defined  in  (4.5). 

Repeating  the  argument  used  to  derive  (4.8)  and  (4.9),  the 
set  in  (4.19)  is  linearly  independent  if 

K<N  (4.21) 

and  it  is  linearly  dependent  otherwise.  It  follows  that  if  K' 
is  the  smallest  value  of  K  for  which  the  Gram  determinant  of 
the  set  in  (4.19)  is  zero,  the  number  of  system  poles  N  is 
given  by 

N  =  K’  -  1  .  (4.22) 

Having  determined  the  number  of  system  poles  by  increasing  K 
until  K'  is  found,  the  system  poles  can  be  identified  by  con¬ 
sidering  the  set 


(4.23) 

-nb.T  -nbjT 

{(Xy1(nT)-y2(nT)) ,•••  ,(XyN(nT)-\yN+1(nT)),  e  ,***,e  }  . 
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Repeating  the  steps  leading  to  (4.16),  the  system  poles  can  be 
found  from  the  roots  of  the  polynomial  in  \ 


N 


l 

i=0 


,  N-i 


^GN+I  +  1^  (  i  +  1,  i+1) )  /  =  0 


(4.24) 


where,  fG,UTil](,,  .  ...  is  the  (i  +  l,i+l)th  cofactor  of  the 

L  N+I  +  l J ( 1  + 1 ,  1+ 1 ) 

Gram  determinant  of  the  set 


-nb-iT  -nb.T 

(y^(nT) , • • •  ,yN+1(nT),e  ,***,e  }  .  (4.25) 

Denoting  the  roots  of  the  polynomial  in  (  4.  24  )  by  , \2 ,  •  •  • ,  xN,  the 
system  poles  are  given  by 


— p ^  =  ln\ ^/T ,  1=1, 2, . . . , N  .  (4.26) 

Equation  (4.24)  thus  permits  the  identification  of  the  system 
poles.  The  only  difference  between  this  result  and  procedure  and 
those  presented  in  Section  4.1.1  is  that,  in  the  present  case, 
knowledge  of  the  input  is  used  and  only  the  system  poles  are 
identified.  Note  that  in  the  pencil  of  functions  method,  Ewen 
[1979]  and  Jain  [1980]  also  use  the  knowledge  of  the  input.  This 
is  further  discussed  in  Section  4.3.  Thus,  the  present  method 
can  be  considered  as  a  discrete  version  of  the  pencil  of  func¬ 
tions  technique  which  uses  shifts. 

In  the  following  section,  use  of  the  input  is  applied  se¬ 
lectively  and  an  iterative  identification  technique  is  defined. 
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4.2 


Iterative  Approach 


The  selective  use  of  the  knowledge  of  the  input  in  the 
identification  procedure  can  serve  as  a  check  for  the  accuracy  of 
the  identified  system  poles  and  in  the  determination  of  the  re¬ 
gions  where  the  poles  are  located.  Moreover,  the  iterative  var¬ 
iation  of  the  sampling  interval  together  with  the  use  as  knowns 
of  previously  identified  poles  can  be  the  basis  of  an  iterative 
identification  method. 

This  approach  can  best  be  illustrated  through  a  simple  ex¬ 
ample.  Suppose  that  the  system  transfer  function  contains  two 
real  poles  -p^  and  -P2  with  P2  much  larger  than  p-^.  Ex¬ 
cite  the  system  by  an  input  consisting  of  a  single  exponential 

x(t)  =  e-l3t. 


Then,  using  (4.4),  the  sampled  output  is  given  by 

yl(nT)  =  cle'nplT  *  c2<TnP2T  * 

First,  choose  a  sampling  interval  T  and  an  input  pole 
that 


(4.27) 

-b  such 


b  <<  p^  <<  P2  (4.28) 

T  -1/b. 

It  follows  that 


bT  -1  (4.29) 

PXT  «  1 
p2  <<  1 


and 
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y^(nT)  »de-nbT  .  (4.30) 

Carry  out  the  identification  technique  without  using  the  know¬ 
ledge  of  the  input  and  identify  a  single  pole  using  the  set 

{ yx (nT ) ,  y2(nT)} 

and  (4.16)  with  N+I  =  l.  From  (4.30)  it  follows  that  this  will 
result  in  the  input  pole  being  identified. 

Iteratively  decrease  T  and  repeat  the  step  above.  When 
the  sampling  interval  is  small  enough  such  that 


PjT  “1  (4.31) 

p2T  <<  1 


the  sampled  output  will,  using  (4.27),  be  approximately  equal  to 


y1(nT) 


c  ^e 


-npjT 


+  de 


-nbT 


(4.32) 


Because  of  the  influence  of  the  system  pole  at  -p^  on  the  out¬ 
put,  the  identification  of  a  single  pole  as  above  will,  no  longer 
result  in  the  input  pole  being  identified. 

The  identification  of  a  single  pole  which  is  not  the  input 

pole  thus  indicates  the  presence  of  a  system  pole(s)  -p^  such 
-np^T 

that  e  1  ,  n=l,2,...,  is  not  negligibly  small.  Then,  this 

pole(s)  can  be  identified  using  (4.24)  where,  knowledge  of  the 
input  is  used. 

A 

Suppose  that  this  procedure  resulted  in  p.  being  identi- 

A 

-np,T 

fied.  The  accuracy  of  p^  ,  that  is  how  well  e  models  the 

system  contribution  to  the  output,  can  be  checked  by  using 
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-np,T 

e  as  a  known  and  identifying  the  input  pole.  That  is,  use 

(4.24)  with  N= 1  and  the  set 

-np.T 

Jy1(nT),  y2(nT),  e  }  (4.33) 

to  carry  out  the  identification.  If,  following  this  procedure, 
the  identified  pole  is  close  to  the  input  pole,  p^  can  be  judged 
to  be  an  adequate  estimate  of  p^.  The  degree  of  closeness  which 
is  required  depends  on  the  amount  of  noise  corrupting  the  mea¬ 
surements.  In  the  noiseless  case,  the  identified  input  pole  is 
judged  close  to  the  actual  input  pole  if  it  is  within  10%  of  the 
true  value.  In  the  higher  noise  level  cases  considered  in  the 
following ,. the  acceptance  region  increases  to  40%. 

This  last  step  can  then  be  repeated  for  decreasing  values 
of  T.  The  identification  based  on  the  set  in  (4.33)  will  give  a 
good  approximation  of  the  input  pole  when  (4.32)  is  an  adequate 
representation  of  the  system  output.  When  this  is  not  the  case, 
that  is,  when  T  is  small  enough  so  that 

P2T  " 1 » 

and  the  system  output  is  given  by  (4.27)  the  input  pole  will  not 
be  identified.  Then,  the  knowledge  of  the  sj stem  input  and  of 
p^  can  be  used  to  identify  the  other  poles  which  contribute  to 
the  system  output.  The  identified  poles  can  again  be  checked  by 
using  them  as  knowns  and  identifying  the  input.  The  above  pro¬ 
cedure  can  be  continued  to  identify  more  poles. 

Note  that  the  iterative  identification  technique  was  illus¬ 
trated  above  for  a  system  which  only  had  two  real  widely  sepa¬ 
rated  poles  and  an  input  which  was  a  single  exponential.  Neither 
of  these  restrictions  is  an  inherent  part  of  the  technique  which 
can  be  readily  generalized  for  complex  conjugate  poles,  grouped 
poles,  and  multi-exponential  inputs. 
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The  iterative  discrete  identification  method  described  here 
has  the  following  advantages: 


•  It  permits  us  to  find  the  regions  where  the  sys¬ 
tem  poles  are  located. 

•  It  provides  a  checking  procedure  for  the  accuracy 
of  the  identified  poles. 

•  It  gives  a  mechanism  for  determining  the  number 
of  system  poles  by  finding  which  poles,  when  used 
as  knowns,  result  in  an  identified  input  pole(s) 
which  is  closest  to  the  true  input  pole(s).  As 
discussed  in  Chapter  5,  this  is  especially  help¬ 
ful  when  the  system  output  contains  noise. 

•  By  using  previously  identified  poles  as  knowns, 
it  allows  the  separate  identification  of  the 
poles  of  a  wideband  system  in  different  frequency 
bands . 


The  only  restriction  inherent  in  this  technique  is  that  the 
system  excitation  be  the  sum  of  exponentials. 

Note  that  if  the  regions  where  the  poles  are  located  and 
the  number  of  poles  in  each  region  are  known  a  priori,  the  itera¬ 
tive  procedure  is  not  even  needed.  For  example,  suppose  that  it 
is  known  that  the  system  contains  three  poles  between  f^  and  f2 
and  two  poles  between  fj  and  f^,  f3>>f2*  The  steps  would  then 
be: 


(1)  Input  e 


-tf' 


(2) 

(3) 


Sample  with  sampling  periof  l/f2 


Use  the  knowledge  of  the 
poles.  The  accuracy  of  the 
be  checked  by  using  them  as 
input. 


input  and  identify  three 
three  identified  poles  can 
knowns  and  identifying  the 
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(  4  )  Input  e 


n 


-tf4 

(5)  Sample  with  period  l/f4. 

(6)  Use  the  knowledge  of  the  input  and  of  the  three  iden¬ 
tified  lower  frequency  poles  to  identify  two  poles. 

(7)  Check  for  the  accuracy  of  the  identified  poles. 

Note  that  this  procedure  would  require  the  factoring  of  a 
polynomial  of  degree  3  to  find  the  lower  frequency  poles  and  of  a 
polynomial  of  degree  2  to  find  the  higher  frequency  poles.  This 
may  result  in  large  errors  in  pole  locations  even  though  the 
polynomials  which  are  factored  are  close  in  a  global  sense  to  the 
true  polynomials.  This  question  is  discussed  further  in  Section 
7.2.  Note  also  that  an  iterative  technique  where  the  poles  are 
identified  one  at  a  time  does  not  require  the  factoring  of  poly¬ 
nomials  since  at  each  step  the  polynomial  is  of  degree  1. 

Before  applying  the  discrete  iterative  method  in  Chapter  5 
to  the  identification  of  a  system  in  the  presence  of  noise,  the 
following  section  relates  the  iterative  method  to  the  pencil-of- 
functions  method  described  in  Chapter  3. 


4 . 3  Relation  to  the  Pe ncil-o f-Funct ions  Method 

Aside  from  the  selective  use  of  the  input  and  the  iterative 
nature  of  the  procedure,  the  iterative  discrete  identification 
method  described  in  Sections  4.1  and  4.2  is  very  closely  related 
to  the  penci 1-of-f unctions  method  presented  in  Chapter  3.  The 
principal  difference  between  the  two  methods  lies  in  the  func¬ 
tions  included  in  the  set  which  is  tested  for  linear  indepen¬ 
dence  . 

More  precisely,  in  the  pencil-o f-f unctions  method  the  poles 
are  determined  using  the  polynomial  equation  in  (3.11)  and  the 
Gram  determinant  of  the  set 
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A 


Iyi(t) . yN+l(t)'  x2(t) . ’Wtl> 

whose  elements  are  formed  by  successive  integrations  as  defined 
in  (3.2). 

On  the  other  hand,  in  the  discrete  iterative  approach  the 
determination  of  the  poles  is  based  on  a  similar  polynomial 
equation  (4.16)  but  with  the  set  now  being  defined  as 

{yx(nT),  ... ,  yN+I+1(nT)} 

or  -nb.T  -nbjT 

{y1(nT),  ...,  yN+1(nT),  e  ,  ...,  e  } 

where,  the  elements  are  formed  by  successive  shifts  as  in  (4.5). 

Note  from  (4.5)  that  successive  shifts  preserve  the  dimen¬ 
sionality  of  the  space  of  y^nT).  That  is,  if  y^(nT)  can  be 
expressed  as  the  linear  combination  of  the  basis  functions 

-np^T  ~npNT  “nb^T  -nbj.T 

e  ,  .  . .  ,  e  ,  e  ,...,e  , 


yj^nT),  obtained  through  (k-1)  successive  shift  of  y^(nT), 
can  also  be  expressed  as  a  linear  combination  of  the  same  func¬ 
tions.  This  space  preservation  aspect  represents  a  fundamental 
advantage  of  the  new  approach.  The  pencil-o f-f unctions  method  as 


described  in  Chapter  3  augments  the  space  of 

~Pi  t 

use  of  successive  integrals.  T*  °  x 


If  e 

the  basis  of  the  system  response,  an  integration 


y -i  ( t )  through  the 

“PNt 

,  e  1  represent 
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e 


t 


/ 

0 


-Pi1 


] 


introduces  the  additional  basis  function  1.  Successive  integra¬ 
tions  introduce  the  additional  components  t,t^,  ...»  . 


The  repeated  integral  and  the  shift  methods  of  forming  the 
elements  of  the  set  are,  by  no  means,  unique.  In  general,  the 
elements  of  the  set  can  be  considered  to  be  the  result  of  a 
mapping  from  the  input/output  functions.  Several  continuous  and 
discrete  domain  mappings  have  been  tried.  Sarkar  [Sarkar,  et 
al.,  (1980)]  suggested  using  reverse  time  integration  which  for 
exponential  inputs  does  not  augment  the  space.  Jain  [Jain  and 
Osman  (1979),  uain  (1980)]  successively  used  the  z-domain  opera¬ 
tor  z/(z-b )  to  form  the  members  of  the  set.  This  technique 
augments  the  space  by  introducing  additional  poles  at  z=b. 


Note  that  the  space  dimension  preservation  property  of 
either  successive  shifts  or  successive  reverse  integrations  de¬ 
pends  on  the  exponential  nature  of  the  input.  If  the  input  can¬ 
not  be  expressed  as  the  sum  of  exponentials,  these  operations 
will  augment  the  space.  Since  an  important  part  of  the  discrete 
iterative  procedure  developed  in  Section  4.2  is  the  ability  to 
selectively  use  the  knowledge  of  the  input  and  since  this  possi¬ 
bility  also  depends  on  the  exponential  nature  of  the  input,  we 
adopted  the  successive  shift  formulation  which  does  not  augment 
the  space. 

Note  also  that  the  selective  use  of  the  input  can  be  ap¬ 
plied  in  the  penci l-o f-f unctions  method.  In  order  to  identify 
the  input  pole  it  would  suffice  to  consider  that  the  input  pole 
is  part  of  the  unknown  system  whereas  the  input  is  an  impulse. 
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SECTION  5 

IDENTIFICATION  IN  THE  PRESENCE  OF  NOISE 


5.1  Introduction 

When  the  identification  technique  is  applied  in  practice  to 
a  real  system,  the  measured  discrete  output  will  be  contaminated 
by  noise.  Then,  instead  of  y^(nT),  the  measured  quantity  can, 
in  general,  be  expressed  as 

z-^nT)  =  y1(nT)  +  WjfnT)  (5.1) 

where  w^{nT)  is  an  additive  noise  which  can  include  quantiza¬ 
tion  noise,  measurement  errors  and  system  thermal  noise. 

Since  the  linear  independence  relations  permitting  the 
identification  of  the  system  poles  will  then  no  longer  be  strict¬ 
ly  valid,  difficulties  can  be  expected  in  the  determination  of 
the  number  of  poles  and  the  identification  of  the  poles  them¬ 
selves.  More  precisely,  the  number  of  poles  is  determined  from 
the  value  of  the  Gram  determinant 


Gj^  =  det 


ziV 


Vk 


.zkzr 


zkzk 


(5.2) 


of  the  set  (z^(nT),  ...,  z^tnT)}  where,  the  shifted  versions  of 
the  measured  output  are  again  defined  as 

z.(nT)  =  zi_1( (n+1 )T) ,  2  <  i  <  k  . 
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and 

z • z  .  =  £  z.(nT)  z*(nT)  .  (5.3) 

1  J  n=0  1  -> 

where  L  is  the  number  of  sampling  periods.  Because  of  the 
noise  term  contained  in  z^(nT),  the  inner  products  which  are 
the  elements  of  the  Gram  matrix  will  also  be  corrupted  by 
noise.  Using  (5.1), 


zizj  =  yiyj  +  yiwj  +  wiyj  +  wiwj  •  (5.4) 

The  first  term  in  the  sum  in  (5.4)  is  the  inner  product  which 
would  be  measured  in  the  absence  of  noise  and  on  which  the  iden¬ 
tification  technique,  as  described  in  Chapter  4,  is  based. 

Knowledge  of  the  statistics  of  the  noise  w^(nT)  permits 
the  calculation  of  the  statistics  of  the  inner  products  and  thus 
can  be  used  to  evaluate  the  effects  of  noise  on  the  inner  prod¬ 
ucts.  The  analysis  of  the  effect  of  noise  on  the  Gram  determin¬ 
ant  in  (5.2),  which  serves  to  determine  the  number  of  poles  and 
is  the  basis  for  the  determination  of  the  poles,  is  highly  com¬ 
plicated  because  of  the  nature  of  a  determinant.  Consequently, 
the  evaluation  of  the  identification  procedure  and  the  determina¬ 
tion  of  the  parameters  to  be  used  in  the  identification  is  best 
carried  out  through  simulation. 


5. 2  Simulation  Program  Characteristics 

A  listing  of  the  simulation  program  used  in  the  evaluation 
of  the  discrete  iterative  technique  and  in  the  determination  of 
the  optimum  identification  parameters  is  given  in  Appendix  A. 
The  functions  performed  by  the  program  and  the  options  available 
can  be  summarized  as  follows. 
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5.2.1  System  Output 

The  program  calculates  the  system  output  as 

z^nT)  =  yj1  ^  (nT)  +  yj2){nT)  +  w(nT),  0  <  n  <  P  (5.5) 

where,  the  output  due  to  one  exponential  at  the  input  is,  using 
(4.3)  and  (4.4),  given  by 


(5.6  ) 


,(j> 


(nT)  = 


R  . 
l 


N 

l  b  -p 
i=l  Dj  Pi 


-np  .T 
e  1  + 


N 


R  . 
1 


. A  .  p . -b  . 
1=1  3 


-nb  .T 
.  3 


j=l,2 


and  w(nT)  is  a  real,  zero-mean,  uncorrelated  Gaussian  random 
process. 

The  number  of  system  poles  N,  the  system  poles 
-p^,  ...  -Pfq,  the  system  residues  R^,  ...  RN,  the  input 
poles  -b^,  -b2 ,  the  sampling  interval  T,  the  length  of  the 
interval  of  interest  L  and  the  standard  deviation  of  the  noise 
are  inputs  to  the  program.  The  choice  of  T  and  L  is  made  at 
each  step  of  the  procedure  in  the  manner  described  in  Sections 
5.3.2  an  d  6.1. 


Several  notes  are  in  order: 


•  Since  the  simulated  system  is  linear,  the  output 
in  (5.1)  which  is  the  sum  of  the  outputs  due  to 
single  exponentials  is  the  same  as  that  due  to 
the  sum  of  the  same  two  input  exponentials. 
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•  For  simplicity,  the  input  poles  are  constrained 
to  be  complex  conjugates  of  each  other.  This 
permits  to  specify  a  real  exponential  input  by 
setting  I^b^}  =  Im < b2 >  =  0;  a  sinusoidal  input 
by  setting  Rgib^}  =  Re(b2}  =  0?  a  decaying  sin¬ 
usoidal  input  by  specifying  and  b2 

(b^=b2*)  as  having  non-zero  real  and  imaginary 
parts . 


•  The  noise  sequence  w^(nT)  models  thermal,  mea¬ 
surement  and  quantization  noise.  For  simplicity, 
it  is  assumed  to  be  Gaussian  even  though 
Wi(nTl,  and  especially  its  quantization  noise 
components  can  be  expected  to  be  non-Gaussian  in 
practice.  If  dithering  is  used,  however,  the 
Gaussian  assumption  can  be  expected  to  be 
valid.  Since  the  s ignal-to-noise  ratio  of 
z^JnT)  is  a  function  of  the  length  of  the  inter¬ 
val  over  which  it  is  defined,  the  noise  standard 
deviation  is  defined  relative  to  the  maximum  value 
of  the  signal  component  yj^(nT)  +  y^^nT) 
instead  of  being  defined  in  terms  of  total  signal 
power  over  the  interval.  This  permits  the  com¬ 
parison  of  results  obtained  using  various  inter¬ 
val  lengths. 


5.2.2  Identification  Algorithms 

In  the  case  where  knowledge  of  the  input  is  not  used  and 
the  input  poles  are  identified  together  with  the  system  poles, 
the  program  performs  the  operations  described  in  Section  4.1.1. 
It 


5-4 


forms  by  successive  shifts  the  set 


{z1(nT) , . .  .  ,zk(nT) } ,  0  <  n  <  L,  L  +  k  <  P  (5.7) 

calculates  the  inner  products 
_  L 

z.z.  =  l  z.(nT)z.(nT) 

->  n=0  -1 


computes  the  Gram  determinant 


Gk  -  det 


z,z.  •••  z.z, 
11  Ik 


Lzkzl  “*  2kZk 


finds  the  (i+1,  i+l)th  co-factors  of  the  Gram  deter¬ 
minant  for  i=0, 1 , . . .  ,k-l 

forms  the  polynomial  equation 


Jo 


1/2  .  0 


finds  the  roots  of  the  polynomial 

determines  the  corresponding  poles  as 
An^i^T'  is  1 1 2  ,  . . . ,  k- 1 
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These  calculations  are  performed  for  a  specified  range  of  values 
of  k . 

In  the  case  where  the  input  and/or  previously  identified 
poles  are  used  as  knowns ,  the  set  in  (5.7)  is  modified  to  include 
these  quantities.  If  a^,a2#.../aj  denote  the  known  input 
and/or  system  poles  the  set  is  defined  as 


-na.T  -na  T 

{ z  ^(  nT ),...,  z^(  nT) ,  e  , .  . .  ,  e  }  •  (5.8) 


The  operations  described  above  are  then  carried  out  for  the  set 
in  (5.8).  Using  (4.23),  the  poles  of  interest  are  found  from  the 
roots  of  the  polynomial 


k-1 


l 

i=0 


k-l-i 


( i+1,  i+1) 


0  . 


5. 3  Parameter  Choices 

Before  carrying  out  the  identification  procedure  it  is 
necessary  to  choose  the  excitation,  the  sampling  rate  or  interval 
and  the  interval  over  which  the  inner  products  which  define  the 
Gram  matrix  elements  are  calculated  so  as  to  insure  optimum  per¬ 
formance  with  a  minimum  of  complexity. 


5.3.1  Excitation 

In  the  choice  of  the  excitation  the  following  factors  must 
be  considered: 
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(1)  The  excitation  must  be  easy  to  implement  in  a  labora¬ 
tory. 

(2)  It  must  excite  the  system  poles.  Since  the  response 
due  to  the  system  is  the  sum  of  decaying  real  or  com¬ 
plex  exponentials  (4.3),  the  identification  is  based 
on  the  transient  system  response.  This  transient  re¬ 
sponse  must  be  excited  by  the  input. 

(3)  It  must  be  adjustable  so  that,  as  a  first  step,  it 
will  essentially  only  excite  the  linear  part  of  the 
system. 

(4)  It  must  result  in  an  easily  analyzable  and  identifi¬ 
able  output. 


An  input  consisting  of  a  sum  of  exponentials  satisfies  the 
four  objectives  stated  above.  We  have  selected  the  following 
elementary  inputs:  a  single  real  exponential;  a  sinusoid;  and  a 
decaying  sinusoid.  They  present  the  advantage  of  containing  one 
or  two  input  poles  which  can  be  readily  identified  using  the  dis¬ 
crete  iterative  approach. 

Although  they  satisfy  requirements  (1),  (3),  and  (4),  the 
decaying  and  nondecaying  sinusoidal  inputs  give  rise  to  the  fol¬ 
lowing  effects  with  respect  to  requirement  (2): 

First,  suppose  that  the  system  contains  one  real  pole  -  p^ 
and  is  excited  by  a  decaying  sinusoid 


x ( t )  =  e“( b+3w) ‘  +  e-(b-j«)t 
=  2e  cosut 

Then,  using  (4.3),  the  system  output  is  given  by 
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R-i  R,  -p,t 

y(t>  ’  ( 


(5.10) 


+  R1  -(b+jo»)t  R1  -  (b- ju )  t 

Pj^-b-jco  p^b+jw 


If  b  is  equal  to  p-^,  the  system  output  becomes 

v(t)  =  e"(b+jal,t  +  ^  -(b-jw)t 

1 '  7  -  jw  ja) 


(5.11) 


and  the  output  does  not  contain  a  contribution  due  to  the  system 
pole.  Thus,  if  it  is  desired  to  identify  p^,  care  should  be 
taken  not  to  use  a  decaying  sinusoidal  input  with  decay  factor 
b  which  can  equal  p^.  Rather,  we  should  select  a  value  of  b 
which  is  smaller  than  the  lowest  real  poles  to  be  identified.  On 
the  ether  hand,  suppose  that  pK  is  one  of  several  system  poles, 
that  Pk+i  is  the  next  larger  pole,  and  that  pK  is  the  largest 
pole  that  has_pa^ready  been  identified.  Then,  if  the  system  is 
excited  by  e  K  cosot,  the  pole  at  -pK  will  not  contribute  to 
the  output.  This  will  facilitate  accurate  identification  of 
pK+l  by  decreasing  the  contribution  to  the  output  of  the  known 
pole  closest  to  the  unknown  pole. 

The  second  effect  that  should  be  noted  arises  in  the  case 
where  a  sinusoidal  excitation  is  used.  Suppose  that 


x(t)  =  2coswt 

=  ejut  +  e"ju)t 


(5.12) 
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Then , 


and,  for  simplicity,  that  the  system  has  one  pole 
using  (4.3),  the  system  output  becomes 


y(t)  =  [ 


R, 


3u-p, 


-Dw-p, 


]e 


-Plfc 


2R,p 


1^1 


2  2 
P].  +  W 


cosot . 


(5.13) 


Note  that  the  system  pole,  p^  contributes  a  decaying  exponen¬ 
tial  to  the  output  whereas  the  input  poles  at  contribute  a 
pure  sinusoid.  It  follows  that  the  input  contribution  will  not 
decay  with  time  and  will  eventually  tend  to  be  greater  than  the 
system  pole  contribution,  thus  obscuring  the  pole  contribution. 
Although  we  give  an  example  in  Section  6.2.4  where  a  pure  sinu¬ 
soidal  input  was  used  to  yield  a  successful  identification,  this 
is  not  generally  true.  In  some  cases,  a  steady  sinusoidal  exci¬ 
tation  can  be  impractical  for  detection  of  poles  in  certain  fre¬ 
quency  regions. 


5.3.2  Correlation  Interval  and  Sampling  Rate 

An  indication  of  the  effect  of  the  sampling  rate  and  the 
correlation  interval  can  be  obtained  from  the  following  consider¬ 
ations. 

Using  (5.4),  examine  the  inner  product  over  an  interval 

[0,L] 


_  _  _  _  L  ? 

Z1Z1  ^  ylyl  +  ylWl  +  Wlyl  +  I  wx(nT)  .  (5.14) 


Suppose  tha'  the  noise  w^(nT)  is  a  zero-mean,  stationary,  uncor¬ 
related  Gaussian  random  sequence  with  variance  o^.  Then  the  ex¬ 
pected  value  of  the  inner  product  is  equal  to 
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(5.15) 


E{ z1z1)  =  y1y1  +  ( L+l ) a  2 

where  y^y^  represents  the  signal  component  and  (L+l)a^  is  the 
bias  introduced  by  the  noise.  Similarly,  the  variance  of  z^z^ 
can  be  expressed  as 


Varfz^z^]  =  2 a  y^y^  +  2(L+l)a 


(5.16) 


Note  that  both  the  mean  and  the  variance  of  the  inner  product 
contain  a  noise  component  which  increases  linearly  with  the  in¬ 
terval  length  L.  It  follows  that  it  is  desirable  to  use  a  short 
correlation  interval  over  which  yjyj  includes  a  large  portion 
of  the  total  signal  power  and  the  noise  contribution  is  as  small 
as  possible. 

The  effects  of  the  sampling  interval  or  rate  can  be  ob¬ 
served  by  considering  the  following,  example.  Suppose  that  the 
noisy  system  output  is  given  by 

z^nT)  =  e-nT  +  e-5nT  +  w^nT)  (5.17) 

and  it  is  desired  to  identify  the  two  poles  1  and  5.  Suppose 
that  a  high  sampling  rate  equivalent  to  ten  times  the  highest 
signal  frequency  (5)  contained  in  (5.17)  is  used.  Then, 

T  =  jq  =  .02  (5.18) 

and 

z  x  ( nT )  =  (e_,02)n  +  (e"#1)n  +  w^nT) 

=  ( . 98 ) n  +  ( . 905  )n  +  wx(nT)  (5.19) 

*  y, (nT)  +  w. (nT) 

1  r-io 


k 


where 


y1(nT)  =  ( . 98 )n  +  ( . 905 )n  . 
Consider  the  Gram  determinant 


|G2|  =  det 


ylyl  +  N11 


Ly2yl  +  N21 


yly2  +  N12 


y2y 2  +  N22  J 


(5.20) 


where,  N^j,  i,j=l,2  represent  the  noise  contribution  to  the  Gram 
determinant  elements.  If  the  correlation  interval  is  suffic¬ 
iently  long  so  that  the  y^y  ^  ,  i,j=l,2,  contain  all  the  signal 
power,  for  pj^  =  -1,  p2  =  -5,  and  T  =  .02,  the  Gram  determinant  in 
(5.20)  is  equal  to 


det 


48.46  +  Nu 
46.41  +  N12 


46.41  +  N12 
44.46  +  N22 


(5.21) 


Note  that  the  signal  components  in  the  matrix  elements  in  (5.21) 
are  within  10%  of  each  other.  This  indicates  that  the  determin¬ 
ant  is  close  to  being  ill-conditioned  and  will  be  severely  af¬ 
fected  by  the  noise  components. 

Suppose  that  a  sampling  rate  equivalent  to  the  highest  fre¬ 
quency  contained  in  (5.17)  is  used.  Then 


and 


z1(nT)  =  ( e~ ' 2 ) n  +  (e  X)n  +  w^nT) 
=  ( . 82 ) n  +  ( -  3 7 ) n  +  w1(nT) 


(5.22) 


Then,  for  a  large  correlation  interval,  the  Gram  determinant  in 
(5.20)  is  equal  to 


|G2|  =  det 


7.08  +  Nn  4.64  +  N12 

_4.64  +  N12  3.08  +  N22_ 


(5.23) 


The  relative  spread  of  the  element  signal  components  in  (5.23)  is 
much  greater  than  in  (5.21).  It  can  be  expected  therefore  that 
(5.23)  will  be  better  behaved  than  (5.21)  in  the  presence  of 
noise  and  would  lead  to  a  more  accurate  identification.  Note 
that  identification  is  not  based  on  the  [G2]  matrix  in  (5.20)  but 
would  be  based  on  the  [G^]  matrix  defined  in  (4.10).  The  [G21 
matrix  is  then  one  of  the  co-factors  of  the  [G^]  matrix.  ‘  Con¬ 
clusions  similar  to  the  above  can  be  drawn  regarding  the  [G3]  ma¬ 
trices  resulting  from  using  a  sampling  interval  of  .02  and  .2. 

The  conclusions  as  to  the  desirable  sampling  interval  and 
correlation  interval  can  be  checked  using  the  simulation  results 
presented  in  Table  5.1.  These  results  were  obtained  for  the  out¬ 
put  given  in  (5.12)  without  additive  noise  and  for  various  levels 
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TABLE  5 . 1 

IDENTIFICATION  RESULTS  FOR  TWO  POLE  CASE 


Pi 

P2 

' 

A 

Pi 

P2 

°/ym. 

me 

1 

2 

.998 

2.0005 

0 

1 

2 

COMPLEX 

CONJUGATE 

.01 

1 

2 

COMPLEX 

CONJUGATE 

.005 

1 

2 

1.074 

1.797 

.001 

1 

2 

COMPLEX 

CONJUGATE 

.001 

5 

1 

2 

1.0049 

1.98 

.001 

5 

1 

2 

1.028 

1.94 

.001 

5 

1 

2 

1.11 

1.78 

.005 

5 

1 

2 

COMPLEX 

CONJUGATE 

.005 

I 

1 

2 

1.05 

1.979 

.001 

05 

1 

2 

1.016 

1.98 

0 

05 

1 

2 

1.026 

2.065 

0.001 

05 

I 

2 

COMPLEX 

CONJUGATE 

0.005 

2 

1 

5 

99988 

5.00013 

0 

2 

1  1 

1 

5 

1.0048 

4.968 

0.001 

2 

|  I 

5 

1.085 

4.727 

0.005 

2 

1 

5 

1.33 

4.2 

0.01 

2 

1 

5 

COMPLEX 

CONJUGATE 

0.01 

05 

1 

5 

COMPLEX 

CONJUGATE 

0.01 

CORRELATION  INTERVAL 
SAMPLING  INTERVAL 
SYSTEM  POLES 
ESTIMATED  POLES 

STANDARD  DEVIATION/MAX.  OF  OUTPUT  OF  ADDITIVE 
GAUSSIAN  NOISE 
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of  additive  noise.  The  noise  present  for  o=0  is  due  to  computer 
round-off  errors.  The  best  results  can  be  observed  to  be  ob¬ 
tained  for  a  small  correlation  interval  length  and  a  large  sampl¬ 
ing  interval  (sampling  rate  approximately  equal  to  the  highest 
frequency  contained  in  the  output). 

For  complex  conjugate  poles,  two  cases  present  them¬ 
selves.  If  the  imaginary  part  of  the  pole  is  smaller  than  the 
real  part  the  conclusions  made  above  as  to  the  desi~  le  sampling 
and  correlation  intervals  remain  true.  If  the  iiu„  ^iiary  part  is 
larger  than  the  real  part  the  sampling  rate  must  be  proportional 
to  the  imaginary  part  in  order  to  prevent  aliasing. 


5.3.3  Noise  Reduction 

In  a  practical  application,  the  noise  level  may,  in  some 
cases,  be  too  high  for  the  proper  operation  of  the  identification 
technique.  Then  several  methods  may  be  used  to  diminish  the  ef¬ 
fects  of  noise. 

If  the  primary  component  of  noise  and  signal  output  are  un¬ 
correlated  and  output  noise  is  zero  mean  useful  methods  of  noise 
reduction  include  the  measurement  and  the  averaging  of  responses 
to  periodic  repetitions  of  the  input  and  the  averaging  of  the 
Gram  matrix  elements  obtained  from  repeated  measurements. 

In  cases  where  the  largest  measurement  error  is  due  to 
quantization  error  induced  by  the  A/D  converter,  dithering  of  the 
A/D  input  signal  could  be  especially  helpful.  This  is  accomp¬ 
lished  by  adding  a  pseudorandom  noise  to  the  output  of  the  black 
box  to  cause  a  uniformly  distributed  amplitude  dither  before  the 
A/D  converter.  Digital  samples  at  the  output  of  the  A/D  conver¬ 
ter  can  then  be  taken  repetitively  and  corresponding  data  points 
averaged  to  eliminate  effects  of  the  dither  and  establish  the 
true  value  of  the  signal.  If  quantization  levels  and  timing  re- 
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mained  constant  for  all  the  repeated  responses,  this  method  could 
be  used  to  improve  the  A/D  accuracy.  Alternatively,  the  timing 
can  be  dithered,  within  a  small  range  so  that  samples  are  first 
taken  on  one  set  of  points  on  the  waveform  and  then  on  another. 

In  this  manner,  the  quantization  error  can  be  averaged  out  to  a 
certain  extent. 

The  noise  effects  will  manifest  themselves  during  the  prac¬ 
tical  implementation  of  the  identification  technique.  The  ade¬ 
quacy  of  the  proposed  noise  reduction  techniques  can  only  be 
truly  ascertained  in  practice  where  constraints  as  to  timing  ef¬ 
fects  are  included  and  the  predominant  noise  sources  are  identi¬ 
fied. 


5. 4  Determination  of  the  Number  of  Poles 

As  described  in  Section  4.1.1  the  number  of  system  poles 
can  be  determined  from  the  order  of  the  lowest  Gram  determinant 
to  vanish.  That  is,  if  K'  is  the  lowest  value  of  K  for  which 
Gk  approaches  zero,  the  number  of  poles  is  given  by 

N  =  K'  -  1. 


In  the  case  where  the  system  output  is  not  corrupted  by  noise, 
this  result  is  valid  as  illustrated  by  the  example  presented  in 
Figure  5.1.  The  parameters  are,  in  this  case: 


System  Poles: 

System  Residues: 
Input  Pole: 

Sampling  Interval  T: 
Correlation  Interval 


.033,  .080 
-.06112,  .036 
.12 
5 

L:  10. 
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In  this  example,  there  are  two  system  poles  and  one  input  pole. 
The  system  output,  therefore,  contains  three  exponentials.  From 
Figure  5.1,  the  sudden  drop  in  the  value  of  the  Gram  determinant 
is  observed  at  K ' =4  and  the  number  of  poles  is  therefore  accur¬ 
ately  determined  as  3.  Note  that  because  of  computer  round-off 
errors  the  Gram  determinant  never  equals  exactly  zero. 

When  the  system  output  is  corrupted  by  noise,  however,  the 
determination  of  the  number  of  poles  from  the  value  of  the  Gram 
determinant  for  increasing  values  of  K  becomes  more  problemati¬ 
cal.  This  is  illustrated  in  Figure  5.2  which  was  obtained  for 
the  same  parameters  as  those  used  for  Figure  5.1  but  with  addi¬ 
tive  noise  whose  standard  deviation  ranged  from  . 1%  to  5%  of  the 
maximum  value  of  the  output.  Clearly,  in  the  presence  of  a  sub¬ 
stantial  noise  component  the  determination  of  the  number  of  poles 
is  very  difficult  and  from  the  results  of  Figure  5.2  would  appear 
to  be  impossible  for  noise  levels  higher  than  .1%  of  the  maximum 
output. 

This  suggests  the  advisibility  of  a  different  criterion  for 
determining  N.  In  particular,  we  explored  the  possibility  of 
identifying  the  system  input,  which  is  known,  as  a  check  on 
whether  the  number  of  unknown  system  poles  has  been  properly 
identified. 

We  believe  that  such  an  approach  is  possible  because  when 
the  number  of  poles  is  incorrectly  specified,  the  poles  which  are 
identified  are,  in  general,  very  different  from  the  true  poles 
and  vary  widely. 

As  discussed  in  Chapter  7,  the  identified  linear  transfer 
function  may  be  relatively  close,  in  a  global  sense,  to  the  true 
transfer  function  even  if  the  location  of  the  identified  poles 
are  quite  inaccurate.  For  this  reason,  we  cannot  be  satisfied 
simply  by  examining  the  linear  transfer  function  and  eventually 
the  acceptability  of  results  will  depend  on  the  effect  of  linear 
pole  location  errors  on  the  identification  of  the  nonlinear 
transfer  functions  themselves. 
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Figure  5.1  Determination  of  Number  of  Poles 
in  Absence  of  Noise  for  Three 
Pole  Output  (points  connected  for 
convenience) 
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SECTION  6 

SIMULATION  RESULTS 


6. 1  Introduction 

In  the  identification  of  a  system,  two  cases  may  present 
themselves.  In  the  first  case,  the  regions  where  poles  are  pres¬ 
ent  and  the  number  of  poles  in  each  region  are  known  a  priori. 
In  the  second  case,  this  information  is  not  available  and  it  is 
desired  to  find  the  poles  in  a  region  of  interest,  let  us  say  be¬ 
tween  frequencies  f^  and  f2*  In  the  first  case,  as  discussed  in 
Section  4.2,  the  iterative  aspect  of  the  identification  may  not 
be  needed.  The  identification  technique  then  becomes  very  simi¬ 
lar  to  the  one  proposed  by  Jain  and  Osman  [1980].  The  additional 
feature  of  the  proposed  technique  is  the  possibility  for  checking 
the  accuracy  of  the  identified  poles  by  taking  them  as  knowns  and 
identifying  the  input  poles. 

In  the  case  where  the  regions  where  the  poles  are  located 
and/or  their  number  is  not  known  and  it  is  desired  to  identify 
the  system  between  frequencies  and  f2,  the  steps  of  the  itera¬ 
tive  discrete  identification  procedure,  as  described  in  Sections 

4.1  and  4.2,  are  as  follows  for  a  real  exponential  input: 


1.  Excite  the  system  with 

...  -b  t 
x ( t)  =  e  , 

b  small  (smaller  than  f^  the  smallest  system  pole  ex¬ 
pected  to  be  found).  Using  a  sampling  interval  in  the 
approximate  range  <  T  <  g-  ,  identify  one  pole  us¬ 

ing  (4.16)  with  N  +  I  =  1.  Since  at  least  one  pole, 
the  input  pole,  contributes  to  the  output,  the  Gram 
determinant  is  not  ill-posed  for  N  +  I  =  1. 
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2.  If  the  identified  pole  is  close  to  -b ,  increase  b,  de¬ 
crease  proportionally  T  and  repeat  (1).  Based  on  ex¬ 
perience  with  several  numerical  examples,  in  this  step 
and  in  step  5  below,  the  identified  input  pole  is 
judged  to  be  close  to  the  known  input  pole  if  it  is 
within  10%  of  the  true  value  in  low  noise  cases  and 
40%  in  higher  noise  cases.  Low  noise  conditions  are 
said  to  exist  if  0  <  a/y  <  .005;  high  noise  condi- 

tions  are  said  to  exist  if  .005  <  a/y  <  .05  . 

-'max 

3.  If  the  identified  pole  is  not  close  to  -b,  use  -b  as  a 
known  and  identify  the  pole(s)  contributing  to  the 
system  output  using  (4.24).  This  can  be  done  in  con¬ 
junction  with  the  information  as  to  the  number  of 
poles  given  by  the  value  of  GK  for  increasing  K. 

4.  Check  the  identified  poles  by  using  them  as  knowns  and 
looking  for  the  input  pole  using  (4.24)  with  N=l. 

5.  When  the  identification  is  found  to  be  satisfactory, 
use  the  identified  pole(s)  as  knowns,  increase  b,  de¬ 
crease  proportionally  T  and  again  look  for  the  input. 

6.  Continue  this  procedure  until  all  poles  in  the  region 
of  interest  are  determined. 

As  noted  in  Section  5.1,  a  decaying  sinusoidal  input,  a 
sinusoidal  input  or  an  input  consisting  of  the  sum  of  exponen¬ 
tials  could  also  be  used. 

In  certain  cases,  the  identification  may  not  be  found  to  be 
satisfactory  for  any  of  the  identified  poles.  This  may  arise  for 
three  principal  reasons.  The  noise  may  be  too  large  for  satis¬ 
factory  identification.  In  this  case,  the  noise  would  need  to  be 
diminished  before  proceeding  with  the  identification.  In  prac- 


6-2 


1 


tice,  this  may  be  accomplished  by  sampling  of  repeated  replicas 
of  the  output  and  averaging  the  samples. 

A  second  cause  for  unsatisfactory  identified  poles  may  be 
poorly  chosen  values  of  b  and  T.  This  would  occur,  for  example, 
if  b  is  too  large  (T  «  1/b).  Then,  too  many  poles  may  be  identi¬ 
fied  at  once  or  the  poles  may  be  much  smaller  than  b,  resulting 
in  a  poorly  conditioned  Gram  determinant. 

Finally,  unsatisfactory  results  may  be  obtained  if  complex 
conjugate  poles,  say  -p^+jw^,  -P]_-jw,  with  wi>>Pi  are  being  iden¬ 
tified  with  T  »  1/Pi*  Sampling  at  the  rate  1/T  would  then  result 
in  substantial  aliasing  leading  to  erroneous  results.  These 
could  be  improved  by  decreasing  T. 

Note  that  these  three  cases  may,  just  as  well,  arise  in  the 
non-iterative  penc il-of-f unc tions  method.  Note  also  that  in  the 
non-iterative  pencil-of-f unc tions  method  the  sampling  rates  to  be 
used  are  not  known  unless  the  regions  where  poles  are  located  are 
known  a  priori.  One  sampling  rate  corresponding  to  the  highest 
frequency  of  interest  f2  cannot  be  used  for  wideband  systems  be¬ 
cause  of  the  large  dynamic  range  of  the  possible  poles  between  f^ 
and  f2.  In  all  cases  where  the  number  of  poles  is  not  known  the 
pencil-of-f unc tions  method  involves  iterations.  Only  one  input 
is  used  but  iterations  are  performed  in  order  to  estimate  the 
number  of  poles. 

The  general  question  of  what  represents  a  satisfactory  or 
unsatisfactory  identified  pole  is  closely  tied  to  the  problem  of 
defining  a  quality  criterion  for  the  identification  results. 
This  question  will  be  addressed  in  more  detail  in  Section  7.  In 
the  present  study,  we  concentrated  on  the  location  of  poles  only 
and  not  location  of  zeros  or  the  value  of  residues.  Both  poles 
and  zeros  are  needed  to  judge  the  accuracy  of  transfer  function 
identification.  In  view  of  this,  the  most  obvious  criterion  of 
accuracy  is,  of  necessity,  simply  the  percentage  mean  square 
error  in  the  value  of  the  identified  poles. 
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6 . 2  Simulation  Results 
6.2.1  One  Pole  Case 


In  order  to  evaluate  the  effect  of  noise  on  the  identifica¬ 
tion  of  a  very  elementary  system  consider  the  one  pole  case  where 
the  impulse  response  is  given  by: 

h ( t )  =  e-t. 

At  time  t=0 ,  excite  the  system  with  a  single  real  exponential 
x ( t )  =  e-bt. 

The  system  output  is  then  given  by 


y(t)  =  ^  [e 


-t 


-bt 
“  e  J 


If  the  input  is  used  as  a  known,  one  pole  needs  to  be  identi¬ 
fied.  If  knowledge  of  the  input  is  not  used,  the  problem  is 
equivalent  to  that  of  identifying  two  poles.  The  number  of  poles 
to  be  identified  was  assumed  known  in  both  cases. 

The  results  of  the  simulation  are  presented  in  Table  6.1. 
A  correlation  interval  of  L  =  10  (10  samples)  was  used  in  all 
cases.  The  identified  parameters  are  denoted  as  b  and  p  .  The 
following  conclusions  can  be  drawn  from  these  results. 

Noise  tolerance  is  considerable  when  only  one  pole  is 
being  identified.  Noise  has  much  more  effect  in  the 
two  pole  case.  In  general,  identification  quality  can 
be  expected  to  be  inversely  proportional  to  the  number 
of  poles  being  identified.  The  best  results  are  ob¬ 
tained  when  b  is  small;  that  is,  for  an  input  ap¬ 
proaching  a  step. 
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TABLE  6.1 

IDENTIFICATION  RESULTS  FOR  ONE  POLE  CASE 
IN  THE  PRESENCE  OF  NOISE 
CORRELATION  INTERVAL  L=10 


°/ymax  b  Pi  b  Pi 


.5  .005 

.5  .01 

.  5  .05 

.2  .005 

.2  .01 

.2  .05 

.002  .005 

.002  .01 

.001  .005 

.OOi  .01 


2. 
2. 
2. 
5. 
5. 
5. 
500  . 
500. 
1000. 
1000. 


1.78 
1.36 
.76 
4.73 
4.2 
1  .7 

474.68 

428. 

949.2 

855.2 


1.11 
±  j  -4 
±  j  1.5 
1.085 
1.33 
±  j  2.9 
8.61 
28.16 
16.3 
55.7 


(b)  Knowledge  of  input  not  used,  b  identified 


In  the  identification  of  two  poles,  noise  tends  to 
generate  complex  conjugate  poles  when  the  two  poles 
are  close  together;  it  results  in  relatively  large  er¬ 
rors  when  they  are  far  apart. 


6.2.2  Wideband  Four-Pole  Transistor  Amplifiers  Circuit 

In  this  section  the  discrete  iterative  approach  is  applied 
to  the  wideband  four-pole  transistor  amplifier  circuit  previously 
studied  by  Jain  and  Osman  [1979].  The  schematic  diagram  of  this 
circuit  is  shown  in  Figure  6.1.  The  transfer  function  of  the 
amplifier  is  equal  to 


H(g)  =  _ 8  ( 1 Q7  )  S2  ( S  -  8000  (IQ6)) _ 

(S  +  .033(106  ))(S  +  . 0  80 (  106  ) ) ( S  +  2 5 . 2  (  1 06  ) ) ( S  +  1205. 1(106)) 

It  contains  two  low-frequency  poles,  one  mid-frequency  pole  and 
one  high-frequency  pole.  The  system  frequency  response  is  shown 
in  Figure  6.2. 

The  iterative  identification  procedure  steps  are  as  follows 
for,  initially,  the  case  where  no  noise  is  added  to  the  output. 

Note  that  the  inner  product  correlation  interval,  L,  is  equal  to 
10  in  all  cases.  It  is  supposed  that  the  normalized  frequency 
band  of  interest  is  between  .001  and  1500. 


1.  Excite  the  system  with  an  exponential  with  a  long  time 
constant,  say,  Tc  =  1000, 


x  (nT) 


- ( .001 )nT 
e 
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V, 


0.01 

Uf 

C3  =  5pf 

9!  =  4 

mn 

g3=40mr 

50.0 

pf 

C4  =  0.01  yf 

9  2  =  1 

mft 

g4  =  0.5  mr 

z  =  l  Kft 


Figure  6.1  wide-band  Transistor  Amplifier  Circuit 
(Jain  and  Osman  [1979]) 
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Use  a  sampling  period  of  T=1000,  so  that  T/Tc=l.  Do  not  use 
x(nT)  as  a  known  in  forming  the  Gram  determinant.  Identifying 
one  pole,  we  find, 

b  =  .00099999 

which  is  very  close  to  the  input.  It  can  be  concluded  that  since 
the  system  poles  did  not  perturb  the  input  time  constant  they  are 
all  larger  than  .001. 


2.  Increasing  the  time  constant  and  decreasing  T,  excite  the 
system  with 


x(nT ) 


- ( . 01  )nT 
e 


with  T=100 .  Identify  the  input  pole  using  (4.16)  with 
N  +  I  =  1;  we  get, 

b  =  .02  . 


3.  Since  .02  is  different  from  .01,  it  is  concluded  that  sys¬ 
tem  poles  perturbed  the  identification  of  input  exponential. 
Look  for  the  dominant  system  pole  (the  one  closest  to  .01)  by 
using  the  same  excitation  and  T  but  now  using  the  input  as  a 
known  in  forming  the  Gram  determinant  in  (4.24)  with  N  +  1. 
Find , 


* 

p^  =  .0328  . 
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Test  by  exciting  the  system  with 


x(nT) 


-  ( .  01  )nT 
e  , 


T=100 


and  using  .0328  in  forming  the  Gram  determinant  and  looking  for 
the  input  time  constant.  Find, 

b  =  .00998 


v/hich  is  close  to  the  input  time  constant  which  is,  of  course, 

A 

known.  It  is  concluded  that  p^  =  .0328  is  a  good  estimate  of 
the  dominant  pole,  that  is,  the  pole  closest  to  b. 

5.  The  procedure  is  now  repeated  at  each  step  looking 

successively  for  the  next  higher  frequency  pole.  The  following 
results  are  obtained: 

i.  Input  exponent  b  =  .05,  T=20,  .0328  used  as  known, 

A  A 

find  b  =  .0012.  Since  b  is  not  close  to  b, 

another  pole  is  present  in  the  vicinity  of  .05. 

ii.  Input  exponent  b  =  .05,  T=20,  .0328  and  .05  used  as 

A 

known s.  Find  ~  .0803. 

iii.  Input  exponent  b  =  .05,  T=20,  .0328  and  .0803  used 

A 

as  knowns.  Find  b=  .04998  which  is  close  to  input 

A 

exponent.  The  estimate  =  .0803  is  thus  assumed 

to  be  correct. 

iv.  Input  exponent  b  =  1.  ,  T  =  l.  ,  .0328  and  .0803  used  as 

A 

knowns.  Find  .  9999.  Since,  b  is  close  to  b,  we 
conclude  that  the  system  contains  no  poles  between 
.0803  and  approximately  3. 
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Input  exponent  b  =  106,  T  =  .l,  .0328  and  .0803  used 

a 

as  knowns.  Find  b  =  7.7.  The  input  is  perturbed  by 
a  system  pole. 

vi.  Input  exponent  10.,  T=.l,  .0328,  .0803,  10.  used  as 

* 

knowns.  Find  =  25.12  . 

vii.  Input  exponent  10.,  T=.l,  .0328,  .0803,  25.12  used 

A 

as  knowns.  Find  10.006.  Thus,  p^  is  judged  to  be 
a  good  estimate  of  the  third  pole. 

viii.  Input  exponent  100.,  T=.01,  .0803,  25.12  used  as 

knowns.  Find  99.  9997. 

Hence,  the  system  has  no  poles  between  25.12  and  approxi¬ 
mately  300.  Note  that  only  .0803,  and  not  both  .0803  and  .0328, 
was  used  as  a  known.  Since  (.0803H.01)  »  (.0328)(.  01)  use  of 
both  results  in  an  ill-conditioned  Gram  determinant. 

ix.  Input  exponent  1000.,  T=.001,  .0803  and  25.12  used 

as  knowns.  Find  600. 

x.  Input  exponent  1000.,  T-.001,  .0803,  25.12,  and  1000 

A 

used  as  knowns.  Find  p^  =  1204.66  .  . 

xi.  The  estimate  can  be  refined.  Input  exponent 

1204.66,  T= .001,  .0803,  25.12  and  1204.66  used  as 

A 

knowns.  Find  p^  =  1205.1  . 

xii.  Input  1000.,  At=.001,  use  1205.1,  25.12,  .0803  as 

A 

knowns.  Find  999.87.  Thus,  p^  is  a  good  estimate 
of  the  fourth  pole. 
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The  procedure  was  repeated  for  the  cases  where  the  output 
was  corrupted  by  additive  zero-mean  white  Gaussian  noise  with 
standard  deviation  equal  to  .05  and  .5  of  the  maximum  value  of 
the  noise  free  output.  The  results  appear  in  Table  6.2. 

For  comparison,  the  results  obtained  by  Jain  and 
Osman  [1979]  in  the  absence  of  noise  are  also  included  in  the 
table. 

Note  that  the  identification  quality  remains  good  for  a 
noise  standard  deviation  equal  to  5%  of  the  maximum  output.  The 
average  percentage  error  in  pole  location  is  .25%,  5.6%  and 

18.5%,  respectively,  for  o/ymax  =  0,  .005,  and  .05  where  the  per¬ 
centage  error  of  the  identified  pole  p  with  respect  to  the  true 
pole  p  is  defined  as 


E(p) 


x  100  . 


The  results  obtained  using  the  discrete  iterative  procedure  with 
an  additive  noise  of  °/ymax  =  *005  are  essentially  of  the  same 
quality  as  those  given  by  Jain  and  Osman  [1979]  in  the  absence  of 
noise. 


6.2.3  One  Real,  Two  Complex  Pole  Case 


Consider  a  system  whose  transfer  function  is  given  by 


TABLE  6.2 

IDENTIFICATION  OF  WIDEBAND  TRANSISTOR  AMPLIFIER  CIRCUIT 


RESULTS  s 


0/¥max 

a 

Pi 

P2 

P3 

a 

P4 

0. 

.0328 

.0803 

25.12 

1205.1 

.005 

.0324 

.0801 

25.5 

1169. 

.05 

.0216 

.0713 

21.05 

1053.2 

JAIN  &  OSMAN 

.0 

.034 

.075 

24.9 

1139.5 

TRUE  VALUES:  p1=.033,  p2=.080,  p3=25.2,  p4=1205.1 
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The  system  contains  three  poles:  one  real  pole  at  -.1625  and  a 

complex  conjugate  pair  at  -.0823  ±  j  .306.  Note  that  the  imagi¬ 
nary  part  of  the  complex  poles  is  larger  than  the  real  part.  The 
simulation  results  for  this  case  are  presented  in  Table  6.3  for  a 
noise  standard  deviation  equal  to  0 ,  .5%,  1%  and  5%  of  the  maxi¬ 
mum  noise  free  output.  The  identified  poles  can  be  observed  to 
be  very  close  to  the  true  values  in  the  first  three  cases.  Only 
for  the  highest  noise  level  do  the  identification  poles  really 
differ  from  the  system  poles.  But,  even  in  this  case  where  the 
noise  level  is  high,  the  results  are  of  reasonable  accuracy.  The 
average  percentage  error  in  pole  location  is,  in  this  case,  .26%, 
1.8%,  4.4%,  and  20.6%,  respectively,  for  o/ymax  =  0,  .005,  .01 

and  .0  5. 

In  order  to  illustrate  further  the  performance  of  the  dis¬ 
crete  iterative  approach  and  to  demonstrate  the  use  of  a  decaying 
sinusoidal  input  it  is  instructive  to  examine  the  steps  used  to 
obtain  these  results.  It  is  supposed  in  this  case  that  the 
frequency  band  of  interest  is  between  .01  and  .15. 

The  system  is  excited  by 

-Ibj+jb -JnT  -[b.-jb  JnT 

x  (nT  )  =  e  +  e 

For  o=0,  the  steps  are: 

Input  T=100,  b2=b2=.005.  Looking  for  two  poles,  find 

A  A 

b^=. 0050002,  b2=. 004999.  Therefore,  conclude  that 
no  system  poles  contribute  to  the  output. 

Input  T=25,  b^=b2=»02.  Looking  for  two  poles,  find 

A  A 

b^=. 019785,  b2=. 019856.  Again,  conclude  that  no 
system  poles  contribute  to  the  output. 
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TABLE  6.3 

IDENTIFICATION  RESULTS  FOR  THREE  POLE  CASE 
(ONE  REAL  POLE  AND  COMPLEX  CONJUGATE  PAIR) 


°/ymax 

Pi 

P2 

P3 

0. 

.1613 

.082352 

+  j  .30607 

.082352 

-  j  .30607 

.005 

.1577 

.08106 

+  j  .3094 

.08106 

-  j  .3094 

.01 

.147 

.  0866 

+  j  .31 

.0866 

-  j  .31 

.05 

i _ 

.135 

.037 

+  j  .361 

.037 

-  j  .361 

i 

TRUE  POLES  ARE:  px  =  .1625,  p2  =  .0823  +  j  .306 

P3  -  .0823  -  j  .306 
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poles 


.0859 


n 


Input  T-5,  bi=b^  =  .l.  Looking  for  two  poles,  find 

A  A  “ 

b^=.02,  b2=.13.  Thus  system  poles  contribute  to  the 
system  output.  However,  looking  for  five  poles,  find 
.16,  .1  ±  j.l,  .0823  ±  j  .306  which  contains  the  input 
poles. 

Conclude  that  .  16,  .0823  ±  j  .  306  could  be  the  system 

As  a  check  input  T=2.5,  b1=b2=. 2  and  use  .16, 

.0823  +  j  .306,  .0823  -  j  .306  as  knowns.  Looking  for 
two  poles,  find  b  =  .19958,  b2  =  .19972.  Thus,  con¬ 
clude  that  .0823  ±  j  .  306,  .  16  represent  the  system 

poles. 


Similarly, 

for  a  =.00  5, 

the  steps  are: 

- 

T=100 , 

b^=b2=. 005. 

Find  bx  =  .00501,  b2 

=  .00501  . 

- 

T  =  25 , 

b  j  =b  2  =  .  02. 

Find  b^  =  .0198,  b2  = 

.0199. 

- 

T-5, 

b^=b 2  —  • 1 • 

Find  b^  =  .02,  b2  =  .13 

• 

- 

T  =  5, 

b  1 =b  2  =  •  !• 

Use  .input  .1  +  j  .  1 , 

.1  -  j. 

1  as 

knowns. 

Find 

.17 

for  N=1 ; 

.074  ±  j  .256  for 

N  =  2 ; 

.18, 

±  j  • 

308  for  N=3. 

- 

T-5, 

t>  1  —t) 2  =  •  ^  9 

use  .17  as  known , 

find 

*1  " 

.0656,  b2  = 

.204.  Conclude  that  . 

17  does 

not 

represent  the  system  poles. 

T  =  5,  b^=b2  =  *l»  use  .  074  +  j  .  256  ,  .  074  -  j  .  256  as 
knowns.  Find  b^=  .0828  ,  b2  =  .107.  Conclude  that 
.074  ±  j  .256  may  represent  the  system  poles. 
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T=5 ,  b^=b2=.l»  use  .18,  .0859  +  j  .308,  .0859  -  j  .308 
as  knowns.  Find  b^  =  .1008  ,  b£  =  .106.  Thus  con¬ 
clude  that  .  18,  .0859  ±  j  .308  represent  the  system 
poles. 

This  same  procedure  applies  to  the  case  where  °/ymax=*01* 
The  noise  level  is,  however,  in  this  case,  sufficiently  high  not 
to  permit  to  properly  check  for  the  accuracy  of  the  identified 
poles  using  a  decaying  sinusoidal  input.  Instead,  using  the  real 
input 


the  steps 

are  as 

follows: 

- 

Input 

T=100 , 

b-.  01. 

Find  b  =  .  0074 . 

- 

Input 

T  =  20 , 

b  =  .  05. 

Find  b=.035. 

- 

Input 

T  =  10 , 

b  =  .  1. 

Find  b=.0198. 

- 

Input 

T  =  10 , 

b  =  .  1. 

Use  knowledge  of  input. 

Find  .135 

- 

Check 

inputting  T=10 

,  b  =  .  1  and  using  .135. 

Find  b  =  . 

A 

119.  Thus,  accept  p^= 

.  135. 

— 

Input 

T  =  2 , 

b  =  .  5  and 

use  .135.  Find  .0288. 

- 

Input 

T=2 , 

b  =  .  5  and 

use  .135  and  .5. 
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Find  .037  ±  j  .361. 


1 


Check  by  inputting  T  =  2,  b  =  .5  and  using  .135, 
.037  ±  j  .361.  Find  b  =  .21.  The  identified  b 
indicates  that  the  identified  poles  are  in  error  but 
is  judged  sufficiently  close  to  b  so  that  they  are  ac¬ 
ceptable. 

Note  that  the  case  treated  in  this  section  is  more  diffi¬ 
cult  than  the  four  pole  case  identified  in  Section  6.2.1.  Al¬ 
though  the  system  only  contains  three  poles,  these  are  not  as 
readily  separable  as  in  the  wideband  four  pole  case. 


6.2.4  Narrowband  Four  Pole  Filter 

As  a  final  example,  consider  the  case  where  the  system 
transfer  function  is  given  by: 


H  (S )  = 


(S+l ) ( S+2 ) ( S+3 ) ( S+4 ) 


which  represents  a  narrowband  filter.  The  example  serves  to  de¬ 
monstrate  the  capability  of  the  identification  technique  when  the 
effects  of  the  various  poles  cannot  be  readily  isolated  through 
use  of  different  sampling  intervals. 

The  identification  results  are  presented  in  Table  6.4. 
Their  accuracy  is  acceptable  for  noise  levels  up  to  5%  of  the 
maximum  system  output.  The  average  percentage  error  in  pole  lo¬ 
cation  ranges  from  2.35%  in  the  noiseless  case  to  13%  in  the  case 
where  °/ymax  =  *05.  Two  sets  of  results  are  presented  for 
a/ymax=*01*  The  first  set  was  obtained  using  a  decaying  sinu¬ 
soidal  input;  the  second  using  a  nondecaying  sinusoidal  input. 
The  identification  was  performed  by  identifying  one  system  pole 
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at  a  time.  The  simultaneous  identification  of  all  four  poles 
would  result  in  complex  conjugate,  not  real,  identified  poles. 

The  procedure  used  to  generate  the  results  in  Table  6.4  can 
be  illustrated  by  examining  the  steps  involved  in  obtaining  the 
results  for  °/ymax=*01  using  the  nondecaying  sinusoidal  input 

x(nT)  =  2  cos  nwT 

=  e~jnwT  +  ejnwT  # 

The  frequency  band  of  interest  is  supposed  to  be  between  .01  and 
5. 


Input  T=100,  w=.005.  Identify  two  poles: 

.000036  ±  j  .0049  -  ±  jw. 

Input  T  =  10,  w=.05.  Identify  two  poles. 

.0013  ±  j  .048  *  ±  jw. 

Input  T=5,  w  =  .l.  Identify  two  poles:  .17, 

.019*  ±  jw. 

Input  T=5 ,  w  =  .l  and  use  knowledge  of  the  input  poles 

A 

±  j  . 1.  Identifying  one  pole,  find  p^  =  .59. 

Test  inputting  T=5,  w  =  .l  and  using  .59  as  known. 

Identify  two  poles:  .0038  ±  j  .098  -  ±  jw.  Thus,  ac¬ 
cept  =  .59. 

Input  T=2.5,  w  =  .2,  use  .59  as  known.  Identify  two 

poles:  .025  ±  j  .196  «  ±  jw. 

Input  T=l.  ,  w  =  .5,  use  .59  as  known.  Identify  two 
poles:  .34  ±  j  .37  *  ±  jw. 
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TABLE  6.4 

NARROWBAND  FOUR  POLE  FILTER 


°/ymax 

A 

Pi 

A 

P2 

p3 

P4 

0. 

1.019 

1.93 

3.09 

3.96 

X 

H 

H 

o 

0 

o 

—  -  - 

.696 

2.01 

2.6 

4.18 

.01(2) 

.59 

2.05 

2.6 

3.95 

0.05 

.7036 

..  —  ■ 

1.815 

3.278 

3.897 

(1)  Decaying  sinusoidal  input. 

(2)  Nondecaying  sinusoidal  input. 
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Input  T=l.  ,  u>=.  5,  use  .59  and  input  poles  ±  j  .5  as 

A 

knowns.  Identifying  one  pole,  find  p^  =  2.05. 

Test  inputting  T=^.  ,  w=.5  and  using  .59  and  2.05  as 
knowns.  Identify  two  poles:  .015  ±  j  .498  »  ±  jto. 
Thus,  accept  =  2.05  . 

The  remaining  two  poles  are  found  by  continuing  this 
procedure. 
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SECTION  7 

IDENTIFICATION  PERFORMANCE  ASSESSMENT 


7. 1  Introduction 

In  order  to  assess  the  practicality  of  the  identification 
technique,  it  is  very  important  to  develop  an  accurate  method  of 
measuring  the  quality  of  the  identification  results.  Two  funda¬ 
mental  approaches  to  this  problem  are  possible:  a  measure  of  the 
accuracy  of  the  identified  parameters  (poles  and  zeros/  residues) 
or  a  global  measure  of  the  difference  between  the  impulse  respon¬ 
ses  or  transfer  functions  of  the  identified  and  actual  systems. 

The  choice  of  the  quality  criterion  is,  of  course,  dictated 
in  practice  by  the  purpose  to  which  the  identification  results 
will  be  used.  In  a  control  systems  application,  the  accuracy  of 
the  identified  parameters  may  be  very  important  in  the  design  of 
a  compensator.  In  a  speech  synthesis  system,  on  the  other  hand, 
an  accurate  identification  of  the  impulse  response  of  the  synthe¬ 
sis  filter  may  be  sufficient. 

Note  that  the  parameter  accuracy  approach  is  much  more  re¬ 
strictive  than  the  global  approach.  An  identified  pole  can  only 
be  accurate  if  it  is  close  to  the  true  pole.  The  representation 
of  an  impulse  response  as  a  function  of  poles  and  residues  is  not 
unique,  however.  An  identified  impulse  response  may  be  accurate 
even  though  the  underlying  poles  and  residues  are  quite  far  from 
the  true  values. 

The  purpose  of  the  linear  system  identification  carried  out 
in  this  effort  is  to  serve  as  a  first  step  to  the  characteriza¬ 
tion  of  the  nonlinear  transfer  functions  of  the  system.  A  final 
choice  of  a  quality  criterion  for  this  problem  cannot  therefore 
be  made  until  the  impact  of  this  choice  on  nonlinear  transfer 
function  identification  is  analyzed. 
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The  present  effort  was  limited  to  the  identification  of  the 
poles  of  the  linear  transfer  function.  As  such,  only  the  par¬ 
ameter  accuracy  of  the  identified  poles  could  be  determined.  A 
choice  between  parameter  and  global  accuracy  measures  depends  on 
their  respective  impact  on  the  identification  of  the  nonlinear 
parts  of  the  system.  Therefore,  a  choice  cannot  be  made  at  this 
time  and  only  the  principles  of  the  parameter  and  global  accuracy 
measures  are  discussed  in  the  following. 


7. 2  Parameter  Accuracy  Measure 

A  reasonable  parameter  accuracy  measure  is  the  ratio  of  the 
difference  squared  between  the  true  and  identified  parameter  to 
the  true  parameter  or  the  percentage  error  in  pole  location. 
Thus,  if  p  and  p  represent,  respectively,  the  true  and  identi- 

A 

fied  values  of  a  pole,  the  parameter  accuracy  of  p  can  be  de¬ 
fined  as: 


p (p>  ■  . 

a  ipi2 

and  the  percentage  error  as: 


A 

p(p)  =  x  100 


(7.1) 


(7.2) 


Note  that  in  the  preceding,  the  error  measure  in  (7.2)  was 
used.  Such  an  accuracy  measure  is  easy  to  implement  during  the 
simulation  testing  of  the  identification  procedure.  It  has  two 
major  disadvantages,  however.  If  the  order  of  the  system  is 
erroneously  determined,  it  does  not  readily  permit  a  measure  of 
the  effect  of  extra  or  missing  poles  and  residues.  Secondly,  in 
the  practical  application  of  the  identification  technique  the 
true  parameter  values  are  unknown.  Consequently,  such  an  ac¬ 
curacy  measure  cannot  be  applied. 
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7. 3  Global  Accuracy  Measure 

The  standard  global  accuracy  figure  measures  the  distance 
between  the  impulse  responses  of  the  true  and  identified  impulse 
responses  normalized  with  respect  to  the  total  energy  in  the  sys¬ 
tem  impulse  response: 


G  ( h  ( t )) 

d 


/ 


h  ( t )  -  h  { t ) 


2dt/  / 


h  { t ) 


(7.3) 


Equivalently,  the  accuracy  figure  in  (7.3)  can  be  expressed  in 
the  frequency  domain  as 


Ga(h(t))  =  Ga(H(f))  =  J 


H ( f )  -  H ( f )  2df/  / | H ( f ) |  a f . 


(7.4) 


Such  a  global  accuracy  figure  is  only  suitable  during  the  simula¬ 
tion  evaluation  of  the  identification  technique.  In  practice  the 
true  impulse  response  or  frequency  response  are  not  known.  The 
true  system  response  can  be  measured  but  not  to  an  impulse  since 
a  perfect  impulse  cannot  be  generated  in  practice.  Measurement 
of  the  true  frequency  response  would  not  be  practical. 

The  global  accuracy  measure  in  (7.3)  can,  however,  be  read¬ 
ily  modified  to  account  for  a  non-impuls ive  input  used  to  measure 
the  true  system  output.  This  is  accomplished  by  weighting  the 
performance  measure  by  the  test  input  g(t).  A  practically  im~ 
plementable  test  input  which  is  as  close  as  possible  to  an  im¬ 
pulse  should  be  chosen.  Note  that  the  test  input  should  be  dif¬ 
ferent  from  the  input  used  to  identify  the  system.  Consequently, 
the  weighed  performance  measure  is  defined  as: 
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Pa(h(t)) 


co 

J  [  [h(t)-h(t )]  *  x(t)]2dt 


CD 

/[ h ( t )  *  x(t)]2dt 


—00 


/  |H(f)  -  H(f  )  | 2  | X (f  )  | 2df 

CD 

a> 

/  | H ( f ) | 2  j  X( f  )j2df 

— ao 


(7.5) 


Equation  (7.5)  can  easily  be  modified  to  measure  the  performance 
of  a  sampled  sygtem: 

l  [  ( h (n )-h (n  ))  *  x(n)j2 

Pg( h(n))  . 


n=-® 


l  [  h  (n )  *  x  (n  )  ] 


/ 

-u 


H(e3“)-H(ei“) j2  |x(ejw)|2dw 


(7.6) 


J  H(e^w)  2  X(e^w)  2dw 

-it 


where  X(n)  is  the  sampled  version  of  x(t)  and  w  is  the 
frequency  normalized  with  respect  to  the  sampling  rate. 

As  discussed  in  Section  7.1,  a  definitive  choice  of  a  per¬ 
formance  measure  must  follow  an  analysis  of  the  effects  of  linear 
system  identification  errors  measured  in  a  parameter  or  global 
sense  on  the  identification  of  the  nonlinear  parts  of  the  sys¬ 
tem.  It  should  be  noted  that  if  a  good  global  performance  mea¬ 
sure  is  sufficient  to  insure  adequate  nonlinear  transfer  function 
identification,  a  question  as  to  the  appropriate  linear  system 
identification  technique  arises  since  a  good  global  performance 
measure  does  not  necessarily  require  that  the  identified  poles  be 
identified  with  an  equally  good  rms  accuracy.  Then,  it  may  be 
sufficient  to  identify  a  system  frequency  response  which  is  close 
to  the  true  frequency  response  and,  from  it,  derive  the  system 
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poles  and  zeros.  The  linear  system  frequency  response  can,  for 
example,  be  estimated  from  steady-state  sinusoidal  measure¬ 
ments.  Such  a  technique,  although  potentially  requiring  more 
measurements,  should  be  more  resistant  to  noise  than  the  tech¬ 
niques  discussed  in  this  report  which  first  estimate  the  system 
poles  and  zeros. 

To  illustrate  the  difference  between  the  global  and  param¬ 
eter  accuracies  and  the  relation  of  these  measures  to  the  identi¬ 
fication  technique,  consider  the  following  example.  Suppose  that 
it  is  known  that  the  transfer  function  of  the  system  is  given  by 


H(s)  = 


PfsT 


(7.7  ) 


that  is,  that  the  system  has  no  finite  zeros.  It  is  desired  to 
identify  H(s).  Suppose  that  the  identification  is  based  on  a 
non-iterative  technique  and  that  the  poles  of  H(s),  that  is,  the 
roots  of  P(s),  are  identified  as  the  roots  of  the  polynomial  of 
the  form  (3.11): 


X  ^G2N+lU+l,  i+1^  /  0 


(7.8  ) 


This  implies  that  P(s)  is  identified  as 


P(S)  =  i=0  S  ^G2N+l^i  +  l,  i+J 


(7.9) 


The  relation  between  the  global  and  parameter  accuracy  measures 
is  now  apparent.  The  global  accuracy  depends  on  the  difference 
between  P( jw)  and  the  polynomial  in  (7.8)  with  \=jw.  The  param¬ 
eter  accuracy  depends  on  the  distance  between  the  roots  of  P(s), 


and  those  of  the  polynomial  in  (7.8).  Note  that  the  two  measures 
are  equivalent  when  one  pole  is  identified  at  a  time,  that  is 
when  N=l. 

In  order  to  ensure  a  good  parameter  accuracy  in  the  general 
case  not  only  must  the  polynomial  in  (7.9)  be  accurate  in  a  glo¬ 
bal  sense  but  the  coefficients  of  the  polynomial  (the  Gram  deter¬ 
minants)  must  also  be  accurate.  Note  that  a  good  global  accuracy 
can  be  obtained  even  if  the  estimate  of  the  number  of  poles  is 
wrong  especially  if  the  estimated  number  is  larger  than  N.  Note 
also  that  the  estimation  of  P( jw)  is  essentially  spectral  estima¬ 
tion.  It  follows  that  an  accurate  identification  in  the  global 
sense  is  much  easier  to  perform  than  an  accurate  identification 
in  a  parameter  sense.  However,  the  adequacy  of  the  global  meas¬ 
ure  in  ensuring  an  accurate  nonlinear  transfer  function  identifi¬ 
cation  is  difficult  to  quantify. 


SECTION  8 
CONCLUSIONS 


The  starting  point  for  our  effort  was  "the  pencil-of-f unc¬ 
tions"  transfer  function  identification  method  originally  pro¬ 
posed  by  Jain  [1974]  for  the  identification  of  linear  systems  and 
later  applied  by  Ewen  [1975,  1979]  also  to  the  identification  of 
nonlinear  transfer  functions  of  a  circuit.  Because  of  the  limi¬ 
tations  and  difficulties  encountered  by  Ewen  [1979]  in  the  pract¬ 
ical  application  of  this  method,  the  objective  of  the  present 
study  was  to  devise  a  modified  approach  to  system  identification 
that  would  lend  itself  more  readily  to  actual  laboratory  measure¬ 
ments. 

As  a  result  of  this  investigation,  the  discrete,  pencil-of- 
functions  identification  (DI)  method  described  in  Chapters  4  and 
5  was  formulated.  Its  main  characteristics  are: 


Completely  discrete  formulation  directly  applicable  to 
sampled  data. 

Selective  use  of  the  knowledge  of  the  input  poles  to 
test  for  regions  where  poles  are  present  and  to  deter¬ 
mine  the  reliability  of  the  identified  poles. 

Iterative  search  for  poles  starting  from  the  lowest 
frequencies  and  progressively  moving  up  to  higher  fre¬ 
quency  poles  with  each  iteration.  Adjustment  of  the 
sampling  interval  T  to  match  the  pole  being  sought. 

Use  of  previously  identified  poles  as  knowns  in  the 
set  defining  the  Gram  determinant  during  the  identifi¬ 
cation  of  additional  poles. 

Use  of  time  shifts  to  create  new  members  of  the  equiv¬ 
alent  to  Jain’s  pencil-of-functions  set. 
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The  advantages  of  the  discrete  iterative  approach  can  be 
summarized  as  follows. 


Repeated  integrations  are  not  required  since  time 
shifts  replace  integration. 

The  accuracy  of  an  identified  input  pole  when  it  is 
far  from  system  poles  gives  a  measure  of  the  best  ac¬ 
curacy  that  can  be  expected. 

Regions  where  system  poles  are  present  can  be  deter¬ 
mined  iteratively. 

The  number  of  poles  to  be  identified  is  determined  by 
stopping  at  the  highest  significant  pole  as  opposed  to 
iterative  predetermination  with  a  Gram  determinant  or 
spectral  analysis. 

Re-identification  of  the  input  poles  while  using  the 
identified  poles  as  knowns  gives  a  measure  of  the  ac¬ 
curacy  of  the  identified  poles. 

The  use  of  identified  poles  as  knowns  is  especially 
useful  for  wideband  systems. 

Each  stage  of  iteration  .requires  the  measurement  of 
approximately  only  twenty  output  samples. 

The  highest  required  sampling  frequency  is  equal  to 
the  Nyquist  rate. 


The  discrete  iterative  identification  (DI)  technique  was 
tested  by  simulation  for  the  case  where  the  output  signal  is  cor¬ 
rupted  by  additive  noise.  The  results  of  these  tests  were  pre¬ 
sented  in  Chapter  6.  The  identification  was  successful  for  both 
narrowband  systems  containing  up  to  four  closely  grouped  poles 
and  for  a  wideband  system  when  output  noise  standard  deviation 
levels  did  not  exceed  5%  of  the  noise  free  maximum  system  out¬ 
put.  Additive  output  noise  level  can  generally  be  decreased  to 
the  desired  level  by  averaging  out  sample  values  in  repeated  ex- 
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periments,  but  an  individual  sample  has  to  be  matched  by  the  num¬ 
ber  of  bits  in  the  A/D  quantizer. 

Based  on  these  results,  it  is  concluded  that  the  DI  method 
should  be  capable  in  practice  of  identifying  systems  containing 
separate  groups  of  four  or  fewer  even  closely  packed  poles  in  the 
presence  of  additive  output  noise  whose  standard  deviation  is  of 
the  order  of  1%  of  the  maximum  system  signal  output.  If  the  out¬ 
put  noise  is  mainly  due  to  quantization,  this  should  permit  use 
of  9  or  10  bit  quantizers.  As  noted  above,  the  hiqhest  required 
sampling  frequency  is  equal  to  the  Nyquist  rate.  This  should 
permit  a  relatively  straightforward  implementation  provided  we 
are  not  dealing  with  a  very  high  Nyquist  rate. 

These  results  compare  very  favorably  with  those  obtained  by 
Ewen  [1979]  who  was  able  to  apply  his  method  to  identification  of 
only  a  two  pole  circuit  and  who  concluded  that  for  his  method,  a 
two  pole  system  required  the  sampling  rate  of  4  to  10  times  the 
highest  frequency  of  the  passband  of  the  system;  that  is,  two  to 
five  times  the  Nyquist  rate,  and  an  A/D  converter  with  a  resolu¬ 
tion  of  not  less  than  16  bits.  Furthermore,  the  identification 
of  a  circuit  containing  a  larger  number  of  poles  would  have  re¬ 
quired  an  A/D  converter  with  a  greater  resolution.  Because  16 
bit  (or  higher)  A/D  converters  cannot  be  clocked  at  high  rates, 
Ewen's  method  was  correspondingly  very  limited  as  to  the  band¬ 
width  and  the  order  (number  of  poles)  of  the  system  which  could 
be  identified.  Moreover,  the  method  had  no  special  provisions 
for  the  identification  of  wideband  systems  to  take  advantage  of 
the  large  spread  in  pole  locations  even  in  the  cases  where  the 
total  number  of  poles  was  relatively  small. 

The  discrete  iterative  technique  described  in  this  report 
has  the  following  comparative  advantages.  Thi  DI  method  is  more 
noise  tolerant  and  this  permits  the  use  of  an  A/D  converter  with 
less  resolution.  Since  the  availability  of  high  speed  10  bit  A/D 
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converters  (output,  noise  less  than  5%)  is  much  greater  than  that 
of  16  bit  converters,  the  DI  method  can  be  employed  at  higher 
sampling  rates.  This  characteristic  added  to  the  fact  that  the 
present  technique  requires  a  sampling  rate  equal  to  the  Nyquist 
rate  implies  that  the  DI  method  will  permit  identification  of 
higher  frequency  circuits.  For  example,  if  the  fastest  commer¬ 
cially  available  16  bit  converter  were  of  rate  125  kHz,  the  maxi¬ 
mum  bandwidth  of  a  system  which  could  be  identified  using  Ewen ' s 
method  would  be  limited  to  125/10  or  125/4,  i.e.,  10  to  30  kHz. 
If,  on  the  other  hand,  a  commercially  available  10  MHz  10  bit  A/D 
converter  were  used  with  the  discrete  iterative  approach,  the 
bandwidth  of  the  system  could  be  as  high  as  10  MHz. 

The  advantage  just  cited  appears  to  offset  the  primary  dis¬ 
advantage  of  the  DI  technique  which  is  that  because  of  its  itera¬ 
tive  nature,  more  input  cases  are  needed.  The  sampling  rate  is, 
however,  no  greater  in  each  case  c.ian  the  Nyquist  rate.  More¬ 
over,  the  iterative  nature  of  the  method  and  the  selective  iden¬ 
tification  of  the  input  poles  as  system  poles  permit  the  identi¬ 
fication  of  wideband  systems  without  the  requirement  that  the  re¬ 
gions  where  the  poles  are  located  and  the  number  of  poles  be 
known  a  priori.  In  defense  of  iteration,  it  should  also  be 
pointed  out  that  the  first  stage  of  the  Ewen  procedure,  which  is 
the  determination  of  the  number  of  poles,  is  itself  also  an  iter¬ 
ative  procedure.  Furthermore,  there  is  little  guidance  to  the 
experimenter  what  sampling  rate  should  be  used  in  that  first 
stage. 

Note,  finally,  that  the  system  identification  approach  dis¬ 
cussed  in  this  report  is  not  being  presented  as  the  ultimate  sol¬ 
ution  for  all  circuits  and  may  still  be  impractical  for  the  case 
of  the  most  complex,  widest-band  circuit.  The  DI  method  is,  how¬ 
ever,  another  step  forward,  hopefully  significant,  for  identifi¬ 
cation  of  circuits  with  relatively  few  poles  and  hopefully  con¬ 
tributes  some  enlightment  for  the  analysis  even  of  circuits  with 
many  poles. 
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APPENDIX  A 
PROGRAM  LISTING 


A0-AU2  593  SI6NATR0N  INC  LEXINGTON  MA  F/6  17/2 

NONLINEAR  SYSTEM  IDENTIFICATION  TECHNIOUE  VALIDATION. <U) 

JAN  82  M  RUOKOf  J  J  BUSS6AN6  F 3060 2-80-C-0 104 

UNCLASSIFIED  A276-4  RADC-TR-81-391  NL 


\ 

\ 


m 


r ir  ti :-identc.ftn 
c  inc«7  jfjcai  jon  mooKAH  a2?* 

C  LESLIE  BATES  11-17-80 
C 

C  km*  PROGRAM 

1NCLUDC  *  J  lit  NIC  .  C*N' 

C 

1XXVI 1 >-o 
J  XXVi  ?)«0 
C 

CALL  IW-UT 
CALL  NQRM1 

PO  600  MSAVE«MMM.MSTOP 
c»0  300  RAN? TR» J , RANF 

JXXt  )  ) »  J  X  XU  f  fi ANPTR*  ?“ 1  ) 

J  XX  (  ?)■]  XXV(RAWPT*»?  ) 

CALL  NOISE 
mm-mSaVE 

Ml »KM- J 

WRJTC  <**2000)  AN 
2000  FORMAT </JX. 'M« ' . 13) 

CALL  hatrxj 

if  (kat2.ne.o)  call  HATRxr 

CALL  COfCTR 
CALL  FOLCStlER) 

300  CONTINUE 

mm-mSaVE 
hT-MM-l 

jr  (RANr.GT.iJ  call  avcrag 
600  CONTINUE 

C 

CALL  EXIT 
ENI. 

c 

SUBROUTINE  input 
PY7E  FNaMEMO) 

complex  EXJ .CXJ.f ACTRJ *FACTR2.SUK1 «  SUM2.SUM3 
INCLUDE  'IOENTC.CHN' 

C 

TYPC  INPUT  FILENAME:* 

accept  llOO.FNAMC 

JJOO  rORMATMOAJ) 

FwamC<  *0  ) • 0 

OPEN  (UNI  Trj  .name  »r  name  ‘TYPE  -  'OLt*'  ) 

READ  (1*1000}  NN.lT.r.K.MMH.MSTOP*  JSMrT 
JOOO  FORMAT (2016) 

READ  (J.2000)  MU1  .PELTAT  .SIC-MA 
*000  FORMAT (20r 16.7 ) 

MU2-C0NjG(nUJ  > 

REAP  (1*2000)  <  OMEGA ( ])•]•!• NN ) 

RCAI*  (  2  *  2000  >  #M  }  >*!■)  •  nn  ) 

REaP  ()*JOOO)  MAT?.  JOI'EF  .  RSNF  .  JOlian 
IF  (RanF.CT.O)  READ  (1  *  1000)  ( 1XXV( I  ). l»l .RANrt?) 

JF  (JOUAN.GT.O)  REAP  (i.rooo)  (W( I > , l«l . JOhan) 

CLOSE  <UN]T«f| 

1  *  1 0  PE  F  4 JPU AN ♦ 2 
IF  IMMM.LT. I)  MMM-1 
IF  (MST0P.LT.MMM)  MSTOP-HMM 
vRITC  (4*4000)  NN.rf .kk.mmm.mSTOP* ISmft 
4  000  FORMAT<  lX**N-#*73/JX.#r»'»J3/JXf,K»'.j4/JX.  'START  J  Mf-  Mb  *  .  1  J* 
1  /IX.  'FINAL  M«  *  t 13. /I  X.  'SHIFT  PARAMETER* '•  I  3 ) 

WRITE  (4.4300)  ML*1  . mu?.  TELTAT  .  SIGMA 
4300  FORMAT  (IX.  'Mill  •  '  .F16.7.  '  .  '  . F  \  6 . 7 .  ?X  »  '  MU*b  '  ♦  F  l  6 . 7 *  '  .  '.F16.7/ 


A- 1 


i  » t  »•«  3vvw  »  »  uni  vim  «  *  *  #  •  J  •  # 

format  <  j  *  ,  *  o*r.C'A«  * »  2( r  it .  7  •  '  •  *.Fi6.7*2x>. 

(/ft. 6».7<r|*. 

WRITE  (4,6000)  (RM).I-l.NN) 

FORMAT u.7,  '  ,  '.FJ6.7.2X). 

4/ix.?x.?<F]6.?* *.Fi6.7.2x>>) 

IF  (MAT?. CP. 0)  WRITE  <4.6S00> 

|F  ^HATJ.WC.O)  WRITE  (4»6600> 

UR]1E  (4,6700) 

FORnAT (  *  •  '  ♦  ‘SHIR  * > 

FORMAT DO* ) 

FORMAT  ('4*.'  liWpJAStl'  MAC*  ELEMENTS  F OR  C-RAh  OCT') 

jf  « ior*cr. co. o. and. jouan.co. 0)  write  (4.7000 

FORMAT  ( J X . 'REGULAR  PEF  0*  I /O  FUNCTIONS') 

IF  C  lOI'CF.NE.O.  AMli.  JOUAN.  CO.  O)  WRITE  «  4 . 7  J  00  ) 

FORMAT <1 X. 'PCF IWE  r (M-J .N)«X1 <N>  Y<M.N)>X?(N) • ) 

IF  (  1  Ot'EF  .NC  .  0.  ANO.  JOUAN.KC  •  0)  WRITE  (4.7200) 
format  <  J  x.  '  PCF  ]  k’E  Y ( M- J  -  1  «  W  )  ■  X  1  <  N  )  V  (  H-  J ,  N  )  -  X?  <  N  >  '  ) 

IF  (JOUAN. WE. 0)  WRITE  (4.0000)  JOUAN, <W( I ), J« 1 , JOUAN) 
FORMAT  MX,  'DEFINE  Y  (  f*  -  J  4  1  •  *  .M*N)«2(N)  '  /  J  )  i  '  J  -  *  *13/ 

IX. ' INRUI  W«* ,?<F16. /♦'»'. Flf. 7. 2X>  . 
</lX.0XO(F!6.7.','»Fl6.7.rX)>) 

IF  (RANF.EO.O)  WRITE  (4,9000) 

IF  (RANF.NE.O)  WRJTE  (4.9100)  RANF ,  ( J XXW(  1  )  .  I  »1  .RANFf 2 ) 
FORMAT  MX.  'NO  Incut  TO  RaWPOm  NUMBER  GENERATOR') 

FORMAT  <)  X.  12.  '  1  Nf'UT  S  TO  RANPOM  NUMhCR  C*C  nc  Rat  OR  :  '  / 

IX, 10< JX. is,  ' . ' .  IS) > 


K  T  —  KK- 1 

C  CALCULATE  X(N) » YH(N) 

DO  200  N-l .RR 

EX-(-J . )ipELTATtN 
EXl-CX*hU' 

EX2-CXXMU2 
XI (N)«EXR(EX1 1 
X2(N)«EXr (CX2) 

p  URITE  (4,1001)  N.xi (NJ.X24N) 

D1001  FORMAT (1 X. '  X  (  '  *  1  3  * ' ) • ' , E 1  4 . 7, '•  'tCl4.7*2X»CJ4.7. 
SUM1 ■<  0 . #0.  ) 

SUM2-(0. .0.  ) 

SUM3* ( 0 . ,0. ) 

DO  100  J • 1 , nn 

FACTR) »R< J )/(MUl-OMCGA(  I  )  ) 
rACTR2»-R<  I  )/<mu2-0mCGa(  \  )  > 

E  X  1  •  (  -I  •  )»nMEC*A(  I  1  *  PC  L  T  AT  *  N 

SUM  lr  SUM  1  ♦(FACTRHFACTRr>»eXF  (EX)  ) 

SUM2-SUM24FACTR) 

SUM3-SUM34^ACTR2 

100  CONTINUE 

YM(N)-SUM1- (XI (W)»SUM2)“(X2(N)»SUM3) 

WRITE  (4,2001)  N »  YM ( M ) 

2001  FORMAT! IX. 'Y(',J3*  '  1  •  '  » £  1  4 . 7  • ',  #»E14.7) 

200  CONTINUE 


IF  (JOUAN. CO. 0)  GO  TR  700 
C  INZT  ARRAYS  YYS  ANp  LAMOA$  FOR  AWERAGING 
700  PO  770  1-1.13 

LAMPAS* I )-( 0  .  .  0.  ) 
no  7S0  J-1.J3 

YTS( 1 . J)-(0. .0. ) 

750  CONTINUE 

770  CONTINUE 

RETURN 
CNP 


c 


COMPLEX  2. EX 
INCLUDE  'IDCMTC.CMN 


c 

c  rjwf« 


CX-«-l  .  >»W<  l  >»PELTaT»n 
2*-EXP<EX  > 

2-2*10. /2D1 V<  J ) 

RETURN 

CMI» 

SUBROUTINE  NORM) 
COMPLEX  2 

INCLUDE  *  JI'ENTC.CMN* 


xxn  -MRSIREML  CXI  <)  >  >  > 

X XP3« APS < REAL (X2<  1  )  >  ) 

YYl»*APS<RCAL  c  VH<  s  )  )  ) 

1*0  150  N-2.K* 

TYPT«AKS(REAL(YH<N) ) ) 

If  (  YYDT  .C'  .  YYD)  YTP-VYPT 

150  CONTINUE 

DO  200  N-)*PP 

XI (N)-X) CN>M0./XXD1 
X2<N)-X2<N)i)0./XXP2 
YM(N J«YNIN)f  JO. /YYP 

D  WRITE  (40001)  W«XJ  (N).M,Xr<N)«N.YH(N) 

P5001  FORMAT  <  J  X • 'X14'»I3» ')«*#C)4.?#* *«E|4.7»'  X2<  ' . 13. ' >.'.E1«  . 7. 

V  \  '.'»EM.7»*  Y<  13.  '  '«  '  .  El  4.  •».  '  *E14.7> 

200  CONTINUE 

IF  <  JOUAW. EO . 0 )  GO  TO  900 
DO  500  I-l.JOUAN 

2«EXP<  <-i  jtui  I >*DELTaT> 

ZPJV( J ) ■ AHS  <  REAL  <  2 ) > 

DO  300  N-J. KK 

7-EXP<  <-1 . ) • WC  J ) »  DEL  T AT  *N) 

22DT-aBSCREali2>> 

IT  l2ZPT.CT.rMVU  ))  2PIVII)-22PT 
300  CONTINUE 

500  CONTINUE 

900  RETURN 

END 

C 

SUBROUTINE  NOISE 

c 

INCLUDE  MDENTC.CHN' 

C 

DP  100  N*1 .BP 

CAU  GAUSS <  Y-OISE  •  S 1  Gma* €••  I XX  ) 

Y ( N  J  •  YM<  N  >  4  YNOISE 

P  WRITE  IA,?001)  YWOJSE.N. Y«N) 

P2001  FORMAT  <  IX.  'YNOISE  -  ' » F 1 5 . 3X •  'Y('*J3»')-'.EM.7.  '  .  '.E 

100  CONTINUE 
RETURN 
END 

SUBROUTINE  MATRXI 

c 

COMPLEX  TI f YJ.SUM1 f pet 
INCLUDE  '1DCNTC.CHN' 

c 

C  FILL  IN  MATRIX 

JPREaa.HH  4  1  -  1  Of’Cr  1 2  -  JOUAN 
DO  *00  1-1. MM 
DO  300  J-I.MM 


>20 

DlOO) 

P  1 
0  1 
200 


no  2 pc  n*i  .m, 

u-nmi-XltHni 

>P  I  I  .LT  .IPftlAk)  00  tO  00 

IP  <  I0PCP  .CO.O.O*  .  1  .OT  .  1  |.R£Ak4  )  t  00  to  50 
IP  < 1 .CO. IRREaki  YI-xi<n> 

K  D.EO.IMIEMtll  TI-X7<N) 

00  tO  PO 
1 > -1 -MM4 JPUAN 
CALL  GCTJ< Y1 , 1 1 .N) 

YJ-Y<N4 J-241SMF1 ) 

IP  < J.LT . 1  BREAK)  oo  tO  >70 

IP  <  lOI'CP. CO. 0.0ft.  J.01 .  1PRCAK4  |  I  00  tO  >00 
IP  < J.CO.lBftCAK)  YJ-XHN) 

IP  <J.C0.I»ft£AK4>  >  Y  Jm x  2 ( N> 

00  10  >20 
JJ- J-MN4 JOUAM 
CALL  OETZ<YJ. JJ.N) 

sum -sum  t  yi  xconjg< yji 

URTtC  <40001)  1.J.N.Y1.YJ.SUM1 

POftMATMX. 'I-'. 17.'  J-'.IJ.'  K-'.J3. 

Yl* • .P>2.». * • • .F12.4.  '  YJ. ' ,r >2.*. *,• .P> 2.4. 
•  sup.  :-'.ri4.6>/t'>P]«.4) 

CONTINUE 
YY< I . Jl-SUN) 

yy<  j.  i  > -con jg < sum  > 

UALUCS  FOR  AUCPA01NG 

YYSO.Jl-YYSUiJ)  4  YYU.J) 

IP  (l.NE.J)  YYS< J. I )-YYS' J.  I)  4  YY<2.I> 

CONTI NUC 
CONTINUE 
00  sss  l-> .nn 

WRITE  <40111)  <YY<J.  Jl.J-l.MH) 

FORMAT! 1X.4<F14.£.','.F)4.6)> 

CONtlNUC 

CALL  OCtftN<-l ,PCT) 

UNITE  44.1000)  PET 

FORMAT  OX.  ’GRAM  PET  E«m  NAN  t  ‘OX.E)4.7,‘,,.E14.71 

RETURN 

ENp 


SUBROUTINE  NAT  R  X  2 


c 

c  till 


COMPLEX  PC  T 
1NCLUPE  *  1  PENT  C  .  CNN  * 

IN  |.l  AGONAL  U1TM  UNHASCP  ELEMENTS 
ualUE«m.*<S1Gna«»J) 
lENp-MH 

IP  OOPCF.NE.O)  ICNP-MN-1 
PO  100  1«1 . IENP 

YY  <1.1  )-YY< 1  . 1  I -VALUE 

CONTINUE 
*'0  555  !•!  .mm 

unite  <4,1111)  < YYI 1 . J) . J-l ,MMI 

FORMAT  <lX.4(F)4.*.','.ri4.4)> 

CONTINUE 

CALL  PEtRMt-1  .PET ) 

WRITE  <400001  PET 

Format  <  i  x .  *  GRam  pet,  UNR1ASCP  mao  ELEMENTS  •  .E)  4.7.  •.£14.7) 

RETURN 

ENP 


subroutine  copcir 


complex  SUMI.AVC 


-•E14.7) 


C  LOOP  r OP  COr  ACTOR 
su"i-<c. >o. > 

DO  500  !•!•«* 

CALL  r>CTRn<  I  .COPACTIJ ) ) 

ua!’e  <4.10001  1. cor ACT*. 

,000  EOR«Al  <  I*.  'COr ACTOR'  .13.  '*  -EI4.7.  '•  *  •£>«•’> 

sun) -SUAI «CorACi  ( I  ) 

500  CONTINUE 

AL>0»SU«1 

UR1TE  <«.ioooi  auC 

JOOC  EORAaT ( 1 X. ' AUEAAGE  Or  COrAdORS  '•£!*. 7. *•*•£! *•/> 

RE  1 URN 
ENP 
C 

SUBROUTINE  FOLESI1ER) 

c 

INCLUDE  -JPEN^C.C"N• 

C  CALCULATE  COE*  0r  POL  YwOAI al  fROA  COEACTORS 
C  Al1  JUS T  AA  AS  NECESSARY 
aa* aa- I  OPE r - JPU an 

Al-AI-lOPCf-JOUAN 

JO  DO  100  1*0. MT 

COER 1 1  i*COr  act  < 14 1 ) 

lr  ICOEP • 1 > .EO. I  0. .o. ) >  00  70  50 

coep  <  1 1  *  (Sort  <cocr  <  i  )M*c«-n»*i> 

50  CONTINUE 

p  UP1TC  <4.10011  l.COEril) 

mooi  format ux .  'Coee <  •  •  II.  ‘  '  «r  15.*. '  • '  »r is.* ) 

ioo  continue 

CALL  ROOTCP<COEr.COr  .AT  .LAMpA.IERl 
If  UER.NE.O-  GO  TO  300 
PO  JOO  1*1 .AT 

A*  REAL  <  L  AAPA  < 1 1 ) 

P»a1aaG<LAM|.a<  1  1  > 

POLE  1 »  «  ATARI  p/A> I /PELT  AT 
FOLCK*<alOS<a»»HP»»2> )/<?.* pel TATI 
UR1TC  <4,10001  1 *laaPa< 2 1 • FOlER.  POLE) 

,000  TORMAT  OX.  'LAAPAC  '  *  13> '  >■'  »El«.7»  '•  '  .El  4 .7.  / 

,  IX  # '  P0LER*' .E14.7. •  F0LEI--.E14.7l 

LAAPASI  1  l-LAAPASl  1  1  4  LAAPAM1 

JOO  CONTINUE 

GO  TO  POO 

300  WRITE  (4.J0001  1ER 

jcoo  format  ox.- error  on  Rooicr  -.ui 

c  FR0A  NOW  ON  USE  AT  1NSTEAP  Or  AA  PECAUSE  ONLY  A-l  ROOTS 
900  CONTINUE 

RETURN 
ENP 


SUPROUT  I  NC  AL'ERaC- 


COmi-lEX  pet 
INCLUPE  *  1  I'EnTC  .  CAN  ■ 

ATT-AT-10PEF- JOUAN 
DO  100  1*1 .ATT 

LAApAI  I  I.LAApASt  I  l/RANr 

A-REAL ILAAPAI I  1 1 

P-Al AAGILAAPAT  I  I  1 

rOLE 1 *• ATARI  P/A) 1 /PELT  AT 

POL  ER*  <  alOG<  a*  t  j<  •  II 1/  <  3 . »  PELT  AT  1 
WRITE  <4,10001  1 .L AAp* 1 1 > . TOLER. POLE  1 
EORAAT  I  l  X"  •  AWC-  LAAPAC  '  .1  1.  ‘  )•' »EI*' 


.£14.7./ 


vv 


Vw** 


c 

PO  250  J-I.HH 

no  ?oo  j- i .hh 

YY(J  ,  j)-YYS<  I  ,  Jl/RANF 

roo  continue 

230  CONTINUE 

1*0  3J5  I-I.HH 

WRITE  (4<nn)  <YY<  1  .  j)  ,  J.J  ,hh) 

mi  FORMAT  <1X*5(C14.7»  *  •  * .E14.7) 1 

335  CONTINUE 

CALL  PCTRM<-I .PET  I 
WRJTC  < 4,2000)  PET 

2000  ro«HAT<l>. 'AVG  GRAM  DETERMINANT # . I  X. El  4. 7.  ' 

call  cofctr 
call  roles 

RETURN 

EN!i 

SUBROUTINE  PC  T  Rn  O  F  L  AG  *  DE  T  ) 

COMPLEX  PET .PCTT <21 
INCLUDE  MPCNTC.CKN* 

c 

IF  ( JFLAG.LT .0)  GO  TO  *00 
C  HOVE  HAT  R  J  X  FOR  COFACTOR  1  INTO  WORKING  MATRIX 
IY-1 

PO  200  I ■ 1 • MT 

JF  < I Y.EO. 1FLAC)  JY-JY41 
JY-l 

PO  100  J-l tHT 

IF  < JY.CP. IFLAG)  JY.JY41 

UTVI 1 . J)-YV( Jt* JY> 

JY-JY41 

100  CONTINUE 

J Y-IY4 l 

200  CONTJNU" 

1h»mt 
GO  TO  700 

C  DOING  WHOLE  HATRJX— HOWE  IT  ALL 

*00  )  M«r  mm 

PO  *50  1-1. MM 
PO  *30  J-l.HH 

WYY(I  .  J>-YYO  ,  J) 

*30  CONTINUE 

*30  CONTINUE 

700  CONTINUE 

p  write  M.ioon 

P1001  rORHATl IX. 'PETRH  MATRIX' ) 

P  PO  753  1-J.IH 

P  NAME  <4.11111  <  UVV  (  I  .  J)  .  J»1  »  JH  ) 

PI  1  1  1  FORMATS  JX.  4  (M  4.*.  *  f  *  .F14.*»1X)1 
P733  CONTINUE 

CALL  CGECO(UYY. 13  -  Im. 1FVT.RC0N.C0F ) 

CD  TYPE  *. 'REVISED  MATRIX' 

CP  no  73*  IM 

CP  UR  IK  *4.1111)  <VYY(  1  .  J)  ,  J-l,  1M> 

CP73*  CONTINUE 

p  WRITE  <4.11121  RCON 

PI  )  1  2  FORmaTOX.  'RCON"  '  .El  4,9) 

CC  IT  < (1 .4RC0N1.E0. 1 . >  GO  TO  *00 

] F  (ACON.CO.O.)  GO  TO  000 
CALL  CGED1 CWYY.13. 1m.  ] PVT . OCT T . COF . 101 
CY»RCAL<PEM<7>  1 
PET-PETT  <  1  1M10.MCK) 

P  WRITE  <4.IJJ3)  PC T ?  < I  I . PE T T < 2 1 . PE T 

PI  11  3  FORMAT  UP.  'PE1-'.3<El4.*,','.c»4.*/)> 
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>pip  ti !-ipentc.chn 
C  COMMON  FOR  I  DEMI  PROGRAM  AJ7a 
C  INPUT 

INTEGER  NN.PP.XIi,MNH.NSTO».MAT2.  ISHri.lODEr 
INTEGER  PANE, JOUAN. IXXU(?0> 

REAL  DELTA!. SIGMA 

COMPLEX  nut . OMEGA < 101 .R< IO).U(  JO) 

C  CALCULATE!! 

INTEGER  NM,MT.KT.IXX<2).IPVl(JO 
REAL  2P1WU0) 

COMPLEX  MU2 

complex  xi <e: 200) . X2< 0:200  .  Y<o: 200 ) • yh)o: 200) 
COMPLEX  YT< J3.13>.yys<13.13).COFnCT<13>.UYYI13.13> 
COMPLEX  LAMDAt 1 3 ) . LAMDAS(13).C0EP  <0: 12) .COr<l<  > 

C  COMMON 

COMMON  NN.PP.KK. MMM. MST0P.MAT2. ISMPT 

COMMON  IOPEF.RANF, JOUANiIXXU 

COMMON  MM.Mt .XT , IXX. IPUT 

COMMON  DELTAT. SIGMA. 2PIU 

COMMON  MU) .OMEGA. R.U.MU2 

COMMON  XI .X2.T.YR.TT.TTS.C0PACT.UYT 

COMMON  LAMDA.LAMDAS.COCF.COP 
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j  Rome  Air  Development  Center  2 

$  RAVC  plant,  and  exeat tea  Kue.aA.di,  development,  tut  and  2 

{  selected  acquisition  programs  In  support  of  Command,  ContKol  £ 

2  Commant-eattowA  and  Intelligence  (C3I)  activities.  Technical  2 

*  a.nd  engineering  support  within  areas  of  technical  competence  ft 

b  <s  provided  to  ESP  PKogKam  Offices  (POs J  and  otheK  ESP  $ 

^  elements.  The  prtnclpal  technical  mission  areas  ate  \ 

k  communications ,  electromagnetic  guidance  and  contxol,  suK- 

5  v^nce  of  gKound  and  aerospace  objects,  Intelligence  data  £ 

*  collection  and  handling,  Information  system  technology,  £ 

s  ionospheric  propagation,  solid  state  sciences,  microwave  2 
i  physics  and  electronic  reliability,  maintainability  and  ft 
S  compatibility.  A 


